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Abstract 

Verification  and  Validation  of  the  Spalart-Allmaras  Turbulence  Model  for  Strand  Grids 

by 

Oisin  Tong,  Master  of  Science 
Utah  State  University,  2013 

Major  Professor:  Dr.  Aaron  Katz 

Department:  Mechanical  and  Aerospace  Engineering 

The  strand-Cartesian  grid  approach  provides  many  advantages  for  complex  moving- 
body  flow  simulations,  including  fully-automatic  volume  grid  generation,  highly  scalable 
domain  connectivity,  and  high-order  accuracy.  In  this  work  the  Spalart-Allmaras  model  is 
implemented,  verified,  and  validated  for  high  Reynolds  number  turbulent  flows  in  a  strand- 
Cartesian  solver.  Second-order  convergence  is  achieved  using  the  Method  of  Manufactured 
Solutions  implying  correct  implementation  of  the  turbulence  model.  By  using  the  NASA- 
Langley  online  resource,  specific  flow  cases  are  validated  with  two  independent  compressible 
codes:  FUN3D  and  CFL3D.  The  strand  solver  is  validated  with  zero-pressure  gradient  flat 
plate  and  bump-in-channel  cases,  and  shows  excellent  agreement  with  FUN3D  and  CFL3D 
for  various  aspects  of  turbulent  flow,  including:  velocity  profiles,  turbulent  viscosity  profile, 
coefficient  of  surface  pressure,  and  drag.  Methods  of  handling  sharp  corners  with  strand 
grids  through  combinations  of  strand  vector  smoothing,  multiple  strands  emanating  from 
a  single  surface  node,  and  telescoping  Cartesian  refinement  into  corner  regions  of  the  near¬ 
body  grid  are  investigated  for  a  NACA  0012  case.  For  standard  viscous  high-aspect  ratio 
grids,  smoothed  strands  with  telescoping  Cartesian  refinement  provide  the  most  accurate 
results  with  the  least  complexity.  Mesh  discontinuities  associated  with  use  of  multiple 
strands  at  sharp  corners  produce  more  error  than  with  smoothed  strands.  With  both 


6 


iv 

strand  approaches  -  vector  smoothing  and  multiple  strands  -  targeted  Cartesian  refinement 
is  critical  to  capture  features  near  sharp  corners  where  strand  grids  alone  are  too  coarse  to 
capture.  Other  results  show  agreement  with  FUN3D  and  CFL3D.  By  using  strand  vector 
smoothing  and  telescoping  Cartesian  refinement,  a  NACA  4412  trailing  edge  separation  case 
is  validated  with  comparison  against  CFL3D  and  FUN3D.  Velocity  profiles  show  reasonable 
agreement  with  CFL3D;  however  implementing  preconditioning  to  the  solver  in  the  future 
may  increase  the  accuracy  of  the  solution. 

(78  pages) 
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Public  Abstract 


Verification  and  Validation  of  the  Spalart-Allmaras  Turbulence  Model  for  Strand  Grids 


The  strand-Cartesian  grid  approach  is  a  unique  method  of  generating  and  computing 
fluid  dynamic  simulations.  The  strand-Cartesian  approach  provides  highly  desirable  qual¬ 
ities  of  fully-automatic  grid  generation  and  high-order  accuracy.  This  thesis  focuses  on 
the  implementation  of  the  Spalart-Allmaras  turbulence  model  to  the  strand-Cartesian  grid 
framework.  Verification  and  validation  is  required  to  ensure  correct  implementation  of  the 
turbulence  model. 

Mathematical  code  verification  is  used  to  ensure  correct  implementation  of  new  algo¬ 
rithms  within  the  code  framework.  The  Spalart-Allmaras  model  is  verified  with  the  Method 
of  Manufactured  Solutions  (MMS).  MMS  shows  second-order  convergence,  which  implies 
that  the  new  algorithms  are  correctly  implemented. 

Validation  of  the  strand-Cartesian  solver  is  completed  by  simulating  certain  cases  for 
comparison  against  the  results  of  two  independent  compressible  codes;  CFL3D  and  FUN3D. 
The  NASA-Langley  turbulence  resource  provided  the  inputs  and  conditions  required  to  run 
the  cases,  as  well  as  the  case  results  for  these  two  codes.  The  strand  solver  showed  excellent 
agreement  with  both  NASA  resource  codes  for  a  zero-pressure  gradient  flat  plate  and  bump- 
in-channel.  The  treatment  of  the  sharp  corner  on  a  NACA  0012  airfoil  is  investigated, 
resulting  in  an  optimal  external  sharp  corner  configuration  of  strand  vector  smoothing  with 
a  base  Cartesian  grid  and  telescoping  Cartesian  refinement  around  the  trailing  edge.  Results 
from  the  case  agree  well  with  those  from  CFL3D  and  FUN3D.  Additionally,  a  NACA  4412 
airfoil  case  is  examined,  and  shows  good  agreement  with  CFL3D  and  FUN3D,  resulting  in 
validation  for  this  case. 


Oisin  Tong,  Master  of  Science 
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Chapter  1 
Introduction 

One  of  the  remaining  bottlenecks  in  computational  fluid  dynamics  (CFD)  is  viscous 
mesh  generation  for  complex  geometry.  With  current  practices,  meshing  experts  can  spend 
days  or  weeks  generating  high-quality  viscous  meshes  around  complex  multi-body  geome¬ 
try,  such  as  rotorcraft.  With  projected  improvements  in  hardware,  the  percentage  of  time 
devoted  to  mesh  generation  using  current  practices  will  only  increase  relative  to  total  sim¬ 
ulation  time.  This  lack  of  automation  places  a  heavy  burden  on  CFD  practitioners  and 
severely  limits  the  practical  use  of  CFD  for  design  and  analysis  of  complex  systems.  Fur¬ 
thermore,  these  complex  systems  require  ever-increasing  mesh  sizes,  for  which  scalability 
becomes  a  greater  issue.  Automating  viscous  mesh  generation  and  maintaining  computa¬ 
tional  efficiency,  while  preserving  spatial  and  temporal  accuracy  are  currently  among  the 
greatest  research  challenges  in  CFD  today. 

The  strand-Cartesian  approach  has  shown  great  potential  to  alleviate  many  of  these 
difficulties  [4-6].  Strand  and  Cartesian  grids  allow  the  possibility  of  fully  automatic  volume 
grid  generation  while  enhancing  scalability  and  the  potential  for  high-order  accuracy.  Near 
solid  bodies,  the  strand  approach  automatically  creates  a  prismatic  mesh  along  “strands” 
emanating  from  pointing  vectors  determined  from  a  surface  tessellation  in  order  to  resolve 
viscous  boundary  layers  and  other  near-body  effects,  as  shown  in  Figure  1.1(a).  Away  from 
solid  bodies,  adaptive  Cartesian  grids  resolve  vortical  shedding  and  wake  features  with  ef¬ 
ficient  high-order  algorithms,  shown  in  Figure  1.1(b).  Due  to  the  robust  and  automatic 
nature  of  the  strand-Cartesian  grid  generation  process,  the  technique  is  easily  extensible 
to  moving-body  problems  for  which  the  grid  can  readily  be  regenerated  at  each  time  step. 
Strand  and  Cartesian  grids  communicate  through  implicit  overset  interpolation  [7-9] ,  which 
is  greatly  facilitated  by  the  fact  that  the  entire  strand-Cartesian  mesh  system  can  be  stored 
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wall  spacing 
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(a)  strand  grid  components 


(b)  strand-Cartesian  grid  for  TRAM  rotor 


Fig.  1.1:  Strand  grid  elements  and  example  strand-Cartesian  grid  system  for  the  TRAM 
rotor. 


on  each  processor  due  to  its  compact  grid  representation.  This  allows  for  self-satisfying 
domain  connectivity  [4]  and  reduces  the  percentage  of  time  needed  for  inter-grid  communi¬ 
cation  [10]. 

The  strand-Cartesian  method  has  proven  effective  for  both  high  and  low  Reynolds 
number  complex  flow  and  geometry  [1,4]  using  various  established  solvers.  A  unique  strand- 
Cartesian  solver  has  been  developed  for  the  strand-Cartesian  approach.  This  work  develops 
this  solver  further  to  handle  turbulent  high  Reynolds  number  flow  for  both  simple  and 
complex  geometry.  Modeling  turbulence  has  been  described  by  some  as  a  “black  hole”  of 
research,  because  by  its  very  nature,  it  is  chaotic  and  unpredictable,  thus  making  its  predic¬ 
tion  and  simulation  difficult.  Nobel  Laureate  Richard  Feynman  famously  described  turbu¬ 
lence  as  “the  most  important  unsolved  problem  of  classical  physics”  [11].  While  simulating 
turbulence  is  more  often  art  than  science,  using  current  turbulence  modeling  methods,  we 
may  achieve  good  results  that  agree  with  experiment  for  many  cases  of  engineering  interest; 
provided  the  correct  methods  are  chosen  for  the  correct  cases. 

Currently,  the  most  widely  used  methods  for  modeling  turbulence  consist  of  Direct  Nu¬ 
merical  Simulation  (DNS),  Large-Eddy  Simulation  (LES),  and  Reynolds-Averaged  Navier- 
Stokes  (RANS)  or  Farve-Averaged  Navier-Stokes  (FANS)  turbulence  models.  The  advance- 
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ment  of  more  powerful  hardware  has  made  solution  methods  such  as  DNS  and  LES  more 
desirable  as  they  provide  accurate  yet  computationally  heavy  solutions.  Basic  RANS  and 
FANS  turbulence  models  provide  a  way  to  predict  flow  regimes  with  some  accuracy  but  with 
far  lighter  computational  workload.  This  approach  is  consistent  with  the  goals  of  extreme 
automation  and  computational  efficiency  targeted  by  the  strand-Cartesian  solver.  A  quick 
analysis  of  high  Reynolds  number  flow  over  complex  geometry  could  potentially  be  made 
in  a  fraction  of  the  time  current  turbulence  and  meshing  methods  take.  This  type  of  quick 
analysis  could  be  used  as  “first  cut”  over  complex  geometry  before  more  accurate  simulation 
methods  such  as  DNS  or  LES  are  used  to  resolve  the  flow  held. 

This  work  aims  to  implement  and  validate  the  one-equation  Spalart-Allmaras  turbu¬ 
lence  model  for  a  strand-Cartesian  solver.  The  current  Navier-Stokes  equations  of  motion 
will  no  longer  be  valid,  and  thus  the  Farve-Averaged  Navier-Stokes  equations  of  motion, 
coupled  with  the  Spalart-Allmaras  equation  will  now  be  implemented  and  solved  in  both 
the  strand  and  Cartesian  solvers.  A  variety  of  case  studies  are  analyzed  in  order  to  verify 
and  validate  the  SA  model  within  the  strand-Cartesian  solver. 

External  sharp  corners,  such  as  those  found  on  the  trailing  edge  of  an  airfoil,  present 
a  challenge  in  meshing  that  is  unique  to  strand-Cartesian  methods.  As  strands  eminate 
normal  to  the  surface  from  a  single  node,  sharp  corners  will  obviously  lack  detailed  reso¬ 
lution,  which  is  critical  for  turbulent  how  in  order  to  capture  how  separation  and  resolve 
the  wake.  Various  methods  such  as  strand  smoothing,  multi-strands,  and  overset  Cartesian 
wake  capturing  meshes  are  investigated  in  this  thesis. 

1.1  Objectives 

The  objectives  of  this  work  are  to: 

•  Investigate  the  new  aspects  of  the  Spalart-Allmaras  turbulence  model,  and  implement 
it  with  FANS  as  a  coupled  system  of  equations  into  the  strand-Cartesian  solver. 

•  Verify  the  mathematical  accuracy  of  the  Spalart-Allmaras  model  in  the  strand  grid 
system.  The  Method  of  Manufactured  Solutions  shall  be  used  to  facilitate  this. 
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•  Validate  this  model  using  the  validation  cases  provided  by  the  NASA-Langley  turbu¬ 
lence  resource.  The  validation  cases  used  are  the;  zero  pressure  gradient  flat  plate, 
bump-in-channel,  NACA  0012  airfoil,  and  the  NACA  4412  trailing  edge  separation. 

•  To  examine  different  external  sharp  corner  meshing  strategies. 

1.2  Outline  of  Thesis 

This  thesis  will  comprise  of  several  chapters.  Chapter  2  contains  the  Literary  Review 
and  will  provide  any  background  information  that  the  reader  may  need  for  understanding 
this  thesis.  Chapter  3  provides  the  reader  with  the  general  theory  and  equations  used  for  this 
work.  Chapter  4  presents  the  numerical  methods  used  in  this  work.  Any  key  information 
to  the  reader  about  discretization  and  implementation  methods  are  given  here.  Following 
this,  the  results  of  the  verification  and  validation  cases  will  be  presented  and  discussed  in 
Chapter  5.  Chapter  6  will  conclude  the  thesis  and  make  recommendations  for  future  work. 
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Chapter  2 

Literary  Review 

This  chapter  provides  the  necessary  background  information,  basic  knowledge,  and 
theory  for  the  work  performed  in  this  thesis.  A  brief  review  of  previous  work  done  on  strand 
grids,  turbulence,  turbulence  modeling,  and  verification  and  validation  is  given  below. 

2.1  Strand  Grids 

Cartesian  Adaptive  Mesh  Refinement  (AMR)  methods  of  grid  generation  can  be  made 
fully  automatic  and  completely  robust  for  arbitrary  geometric  configurations,  however, 
Cartesian  AMR  methods  cannot  be  used  exclusively  to  satisfy  HPC  requirements  for  au¬ 
tomation  due  to  the  length  scales  associated  with  viscous  flow.  The  reason  being,  that 
Cartesian  methods  require  isotropic  refinement  to  accommodate  arbitrarily  complex  bound¬ 
ary  surfaces.  Fluid/aerodynamic  interactions  with  solid  surfaces  result  in  viscous  boundary 
layers  that  require  spacing  too  fine  to  be  accommodated  isotropically  [4],  It  is  in  these 
viscous  boundary  layers  that  strand  grids  are  employed  to  provide  the  necessary  grid  reso¬ 
lution,  resulting  in  the  strand-Cartesian  approach. 

The  strand-Cartesian  grid  approach  provides  many  advantages  for  complex  moving- 
body  flow  simulations,  including  fully-automatic  volume  grid  generation,  highly  scalable 
domain  connectivity,  and  high-order  accuracy.  Strand  grids  are  created  from  discrete  surface 
geometry  comprised  on  n-sided  poplygons.  Tessellated  geometries  can  be  formed  from  any 
n-sided  polygon,  including  triangles  and  quadrilaterals  in  3D.  A  strand  is  created  for  each 
vertex  of  the  discretized  surface  geometry.  A  unit  normal  vector  is  approximated  at  each 
surface  vertex  by  averaging  the  normals  of  neighboring  polygons.  Each  strand  is  placed 
so  its  origin,  or  root,  coincides  with  the  vertex  and  initially  points  in  the  direction  of  the 
corresponding  surface  normal.  Every  strand  in  the  mesh  shares  a  length  and  user-defined 
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ID  nodal  distribution  along  the  strand.  Each  strand  has  a  clipping  index  to  “shorten” 
the  strand  to  prevent  intersection  with  other  bodies.  This  is  to  say  that  nodes  that  lay 
beyond  the  clipping  index  are  not  valid  physical  nodes  [1,4].  Fig.  1.1(a)  demonstrates  this 
procedure  on  a  simplified  surface  triangle. 

Sharp  corners,  both  internal  and  external,  present  challenges  for  strand  grids.  At 
internal  corners,  strands  tend  to  intersect  and  penetrate  the  surface  geometry.  Strands 
typically  fail  to  have  sufficient  resolution  at  external  corners.  Strand  smoothing  can  help  to 
mitigate  the  challenges  presented  by  corners.  Strand  smoothing  refers  to  a  simple  method 
of  adjusting  the  surface  normals  to  properly  cover  all  regions  of  the  domain.  Fig.  2.1(a) 
shows  a  completely  un-smoothed  strand  mesh  around  a  square  cylinder  where  Fig.  2.1(b) 
illustrates  this  same  mesh  after  being  smoothed. 

An  alternate  approach  to  handling  sharp  corners  has  been  investigated  by  Work  [12] 
et.  al.  Multi-strands  work  by  placing  multiple  stands  on  a  single  node  at  the  point  of  the 
corner,  and  fanning  strands  evenly  around  this  angle.  Fanning  strands  in  an  even  distibution 
around  a  sharp  corner,  while  simple  in  2D,  provides  a  challenge  in  3D.  In  3D,  the  number  of 
strands  required  from  a  single  node,  and  the  small  angles  that  would  separate  the  strands 
in  order  to  provide  the  neccesary  resolution  at  the  corner,  provides  numerous  numerical 
and  meshing  issues  such  as;  small  unusual  shaped  high  aspect  ratio  cells  with  many  faces. 
The  NACA  0012  case  in  Chapter  5  presents  some  of  these  results  for  a  2D  airfoil.  Strand 
smoothing  with  Cartesian  refinement  in  the  corners  was  shown  to  be  the  most  effective 
method  for  resolving  external  corners. 

Strand  smoothing  is  typically  insufficient  for  sharp  internal  corners.  Even  with  smooth¬ 
ing,  strands  often  intersect  creating  a  negative  volume  element.  Strand  clipping  is  used  to 
ensure  that  strand  cells  do  not  intersect  or  penetrate  the  surface  geometry.  An  example  of 
a  smooth  and  clipped  internal  corner  is  in  Fig.  2.2. 

Strand  grid  technology  has  an  advantage  over  other  mesh  technologies  as  it  collapses 
memory  requirements  for  near-body  domain  partitions  from  volume  data  to  surface  data 
only.  Non-surface  grid  related  data  that  may  be  required  during  a  simulation  is  trivially 
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(a)  No  Smoothing  (b)  Significant  Smoothing 


Fig.  2.1:  Sample  strand  meshes  around  a  square  cylinder  to  illustrate  how  smoothing  affects 
strands  [1]. 

computed.  It  is  easy  to  undervalue  the  advantage  of  minimal  memory  requirements  for  grid 
related  data.  After  all,  memory  is  relatively  cheap.  The  advantage  is  realized  in  context 
of  a  design  environment  with  HPC  systems.  Maximum  solution  throughput  is  obtained  by 
allocating  the  minimum  number  of  processors  that  can  accommodate  a  given  problem  in  core 
memory.  Multiple  cases  can  then  be  run  simultaneously,  realizing  perfect  parallelism  [4]. 

2.2  Turbulence 

In  1883  Reynolds  [2]  published  the  results  from  a  series  of  experiments  in  order  to  char¬ 
acterize  the  different  flow  regimes.  The  experiment  consisted  of  dye  in  a  pipe  that  contained 
different  speeds  of  water  flowing  through  it.  In  slow  water  the  dye  stayed  streamline.  As 
the  speed  of  the  water  increased,  the  dye  became  progressively  unstable.  Figure  2.3  shows 
a  diagram  of  these  findings. 

Figure  2.3(a)  shows  the  streamline  laminar  flow,  while  Figure  2.3(b)  shows  the  eventual 
turbulent  flow  that  occurs  once  a  critical  velocity  has  been  reached.  Reynolds  also  observed 
that  if  the  critical  velocity  was  slowly  reached,  there  was  an  emergence  of  “eddies,”  a  key 
feature  to  turbulent  flow,  shown  in  Figure  2.3(c). 
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Fig.  2.2:  Example  of  internal  corner  with  smoothing  and  clipping  [1], 


(a)  Laminar  Flow  (b)  Turbulent  flow 


V 

(c)  Turbulent  Eddies  within  the  flow 

Fig.  2.3:  Observations  from  Reynolds’  Dye  and  Pipe  Experiment  [2]. 


The  end  result  of  these  experiments  was  the  introduction  of  a  dimensionless  quantity; 
the  Reynolds  number  (Re).  It  is  defined  as  the  ratio  of  inertia  forces  to  viscous  forces. 
It  is  used  to  characterize  different  flow  regimes  within  a  similar  fluid,  such  as  laminar  or 
turbulent  flow:  laminar  flow  occurs  at  low  Reynolds  numbers,  where  viscous  forces  are 
dominant,  and  is  characterized  by  smooth,  constant  fluid  motion;  turbulent  flow  occurs  at 
high  Reynolds  numbers  and  is  dominated  by  inertial  forces,  which  tend  to  produce  chaotic 
eddies,  vortices  and  other  flow  instabilities. 

There  is  no  general  theorem  that  relates  Reynolds  number  and  turbulence.  For  each 
type  of  flow,  a  different  Reynolds  number  is  required  for  turbulence.  For  example,  Poiseuille 
flow  with  a  Reynolds  number  value  over  2300  begins  to  become  turbulent.  Turbulence  flow  is 
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generally  interspersed  with  laminar  flow  until  a  larger  Reynolds  number  of  about  4000  [13]. 

In  spite  of  how  common  turbulence  is,  it  is  not  easy  to  define  precisely,  however, 
turbulence  has  certain  characteristics,  that  Kundu  [14]  describes  as: 

•  Randomness:  Turbulent  flow  is  irregular,  unpredictable  and  chaotic. 

•  Nonlinearity :  Turbulent  flows  are  highly  nonlinear.  The  nonlinearity  serves  two  pur¬ 
poses.  First  it  causes  the  relevant  nonlinearity  parameter  to  exceed  critical  value. 
In  unstable  flows  small  perturbations  grow  spontaneously  and  frequently  equilibrate 
as  finite  amplitude  disturbance.  On  further  exceeding  the  stability  criteria,  the  new 
state  can  become  unstable  to  more  complicated  disturbances,  and  the  flow  eventually 
reaches  a  chaotic  state.  Second,  the  nonlinearity  of  a  turbulent  flow  results  in  vortex 
shedding,  a  key  process  by  which  three-dimensional  turbulent  flows  maintain  their 
vorticity. 

•  Diffusivity :  Turbulent  flows  have  a  rapid  rate  of  diffusion  of  heat  and  momentum  due 
to  the  macroscopic  mixing  of  fluid  particles. 

•  Vorticity:  Turbulent  flows  are  characterized  by  high  levels  of  fluctuating  vorticity.  The 
identifiable  structures  in  a  turbulent  flow  are  vaguely  called  eddies.  Flow  visualization 
of  turbulent  flows  shows  various  structures  -  coalescing,  dividing,  stretching,  and  above 
all  spinning.  A  characteristic  feature  of  turbulence  is  the  existence  of  an  enormous 
range  of  eddy  sizes.  The  large  eddies  have  a  size  of  order  roughly  equivalent  to  the 
width  of  the  region  of  the  turbulent  flow;  in  a  boundary  layer  this  is  the  thickness  of 
the  layer.  The  large  eddies  contain  most  of  the  energy.  The  energy  is  handed  down 
from  large  to  small  eddies  by  nonlinear  interactions,  until  it  is  dissipated  by  viscous 
diffusion  in  the  smallest  eddies,  whose  size  is  of  the  order  of  millimeters. 

•  Dissipation:  The  vortex  stretching  mechanism  transfers  energy  and  vorticity  to  in¬ 
creasingly  smaller  scales,  until  the  gradients  become  so  large  that  they  are  smeared  out 
(i.e.  dissipated)  by  viscosity.  Turbulent  flows  therefore  require  a  continuous  supply 
of  energy  to  make  up  for  the  viscous  losses. 
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In  1941  Kolmogorov  [15]  postulated  that  for  very  high  Reynolds  numbers,  the  small 
scale  turbulent  motions  are  statistically  isotropic,  thus  Kolmogorov  length  scales  were  in¬ 
troduced.  Kolmogorov  length  scales  are  the  smallest  scales  in  the  spectrum  that  form  the 
viscous  sub-layer  range.  In  this  range,  the  energy  input  from  nonlinear  interactions  and 
the  energy  drain  from  viscous  dissipation  are  in  exact  balance.  The  small  scales  are  in  high 
frequency,  which  is  why  turbulence  is  locally  isotropic  and  homogeneous. 

Laminar  flow  may  easily  be  modeled  by  simply  solving  the  Navier-Stokes  equations, 
turbulent  flow  provides  a  challenge  however.  Using  ordinary  CFD  techniques  (not  DNS 
or  LES),  it  is  possible  to  find  some  particular  solutions  of  the  Navier-Stokes  equations 
governing  fluid  motion,  however,  all  such  solutions  are  unstable  to  finite  perturbations  at 
large  Reynolds  numbers.  Thus  different  CFD  methods  must  be  used  (DNS  or  LES),  or 
some  type  of  turbulence  modeling  must  be  utilized. 

2.3  Turbulence  Modeling  and  One-Equation  Turbulence  Models 

2.3.1  Turbulence  Modeling 

Direct  numerical  simulation  is  potentially  the  most  accurate  method  of  modeling  tur¬ 
bulence.  DNS  is  a  CFD  simulation  that  solves  the  Navier-Stokes  equations.  This  requires 
that  the  entire  range  of  temporal  and  spatial  scales  of  the  turbulence  to  be  solved.  In  the 
computational  mesh,  all  spatial  scales  must  be  resolved  [16];  from  the  smallest  dissipative 
scales  (Kolmogorov  scales),  to  the  integral  scale  L,  which  is  associated  with  the  motions 
containing  most  of  the  kinetic  energy.  Due  to  range  of  scales  that  need  be  resolved,  DNS 
requires  powerful  hardware  to  simulate. 

Because  of  the  computational  cost,  DNS  is  used  as  a  tool  in  the  fundamental  research 
of  turbulence.  DNS  is  used  to  better  understand  the  physics  of  turbulence  through  the  use 
of  “numerical  experiments.”  Generally  the  information  extracted  is  difficult  or  impossible 
to  obtain  in  a  laboratory.  DNS  may  also  be  used  in  the  development  of  turbulence  models 
for  practical  applications,  such  as  those  used  in  the  subgrid  in  LES  models,  and  methods 
that  solve  the  RANS  equations,  such  as  the  SA  model. 
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A  practical  and  accurate  tool  for  simulating  turbulent  flows  is  large-eddy  simulation 
(LES).  In  turbulent  flow,  the  large  eddies  are  dependent  on  geometry,  while  the  smaller 
scales  are  more  universal  [15].  The  velocity  field  in  turbulent  flow  may  be  separated  into 
resolved  and  subgrid  components.  The  resolved  portion  of  the  field  represents  the  “large 
eddies,”  while  the  subgrid  part  of  the  field  represents  the  “small  eddies,”  whose  effect  on  the 
resolved  field  is  included  through  the  subgrid  scale  model.  This  feature  of  turbulent  flow 
allows  for  large  eddies  to  be  solved  explicitly  in  calculation.  The  small  eddies  are  implicitly 
accounted  for  by  using  a  subgrid-scale  model.  Subgrid-scale  (SGS)  modeling  [17]  is  used  to 
represent  the  effects  of  unresolved  small-scale  fluid  motions  (small  eddies,  swirls,  vortices) 
in  the  equations  governing  the  large-scale  motions  that  are  resolved  in  computer  models. 

Both  DNS  and  LES  provide  the  most  accurate  results  when  simulating  turbulence, 
but  at  the  expense  of  an  extremely  high  computation  workload.  A  computationally  lighter 
alternative  to  DNS  and  LES  is  to  use  a  RANS  turbulence  model.  A  turbulence  model  aims 
to  model  all  the  unsteady  turbulence  eddies.  No  attempt  is  made  to  resolve  the  unsteady 
features  of  any  of  the  turbulent  eddies,  not  even  the  largest  ones.  Instead,  mathematical 
models  are  employed  to  take  into  account  the  enhanced  mixing  and  diffusion  caused  by 
turbulent  eddies.  When  using  a  turbulence  model,  the  Navier-Stokes  equations  are  no 
longer  employed,  instead  the  Farve-Averaged  Navier-Stokes  equations  are  used  (Chapter 
3). 

2.3.2  One-Equation  Turbulence  Models 

One-equation  turbulence  models  solve  one  turbulent  transport  equation  in  conjunction 
with  the  RANS  or  FANS  equations.  One  turbulent  transport  term,  usually  the  turbulent 
kinetic  energy,  is  solved  by  the  turbulent  equation.  The  original  one-equation  model  is 
Prandtl’s  one-equation  model  [18].  The  Spalart-Allmaras  model,  which  is  the  focal  model 
in  this  work,  is  based  on  Prandtl’s  original  one-equation  model.  Other  one-equation  models 
exist,  such  as  the  Rahman-Siikonen-Agarwal  Model  and  the  Baldwin-Barth  model,  however 
these  do  not  share  the  popularity  of  the  SA  model,  especially  for  external  aerodynamics. 

The  standard  Spalart-Allmaras  model  [19]  consists  of  a  single  transport  equation  which 
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calculates  the  kinematic  eddy  viscosity  (D).  This  model  has  been  shown  to  give  good  results 
for  adverse  pressure  gradients  and  boundary  layers.  Thus  this  model  has  become  popular  for 
aerodynamic  and  turbo-machinery  purposes  [20].  The  original  Spalart-Allmaras  equation 
has  the  basic  form  of: 


p^  +  C  =  F  +  P-D  +  T  (2.1) 

Here  we  define  C  as  the  convection  of  D,  F  as  the  diffusion  of  h,  P  as  the  rate  of 
production  of  v,  D  as  the  rate  of  dissipation  of  z>,  and  T  as  another  term  for  the  transport 
of  v  by  turbulent  diffusion.  For  convenience,  we  may  let  P  —  D  +  T  =  S,  and  call  these 
the  source  terms.  Chapter  3  provides  a  more  explicit  explanation  of  the  model  and  it’s 
modifications. 

The  advantage  to  a  one-equation  model  is  simply  that  there  is  only  one  additional 
equation  to  solve.  Multi-equation  models  require  two  or  more  additional  equations  to  be 
solved.  While  multi-equation  models  may  or  may  not  produce  better  results,  there  is  added 
complexity  when  implementing  the  model  into  the  CFD  solver,  such  as  linearizing  the  source 
terms.  Adding  multiple-equation  turbulence  models  to  the  equations  of  motion  adds  extra 
time  required  for  the  solution  to  converge.  Chapter  4  describes  the  methods  required  for 
discretizing  and  implementing  the  turbulence  model. 

2.4  Verification  and  Validation 

2.4.1  Verification 

Whenever  new  numerical  algorithms  are  designed  and  implemented  in  a  CFD  code, 
it  is  necessary  to  perform  detailed  verification  to  ensure  accuracy  and  to  test  the  code’s 
integrity.  Verification  in  this  sense  is  defined  as: 

The  process  of  determining  that  a  model’s  implementation  accurately  represents  the  devel¬ 
oper’s  conceptual  description  of  the  model  and  the  solution  to  the  model.  [21] 
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Roache  [22],  Roy  [23],  and  Veluri  [24]  showed  how  the  methods  of  manufactured  so¬ 
lutions  (MMS)  can  be  used  to  verify  the  desired  order  of  accuracy  of  the  scheme  with  the 
actual  order  of  accuracy  the  code  produces.  They  demonstrated  that  verification  using 
MMS  is  capable  of  identifying  malformed  algorithms  and  general  coding  “bugs.” 

Many  researchers  have  extended  MMS  to  verify  complex  problems.  Diskin  et  al.  [25-27] 
have  employed  MMS  extensively  to  compare  and  test  different  second  order  schemes  as  well 
as  study  the  effect  irregular  grids  have  on  accuracy.  Many  others  have  developed  MMS 
methodologies  for  the  interior  solution  only  [28-33] ,  while  Folkner  and  Katz  [34]  investigated 
this  methodology  on  boundary  conditions. 

Roache  [22]  describes  a  general  MMS  provides  a  general  procedure  for  working  with  such 
analytical  solutions.  The  procedure  is  very  simple.  A  continuum  solution  is  constructed, 
which  in  general  will  not  satisfy  the  governing  equations.  An  appropriate  source  term  can 
be  determined  to  cancel  any  imbalance  in  the  PDEs  caused  by  the  choice  of  the  continuum 
solution.  The  solution  also  defines  the  boundary  conditions  in  all  forms,  be  they  Dirichlet, 
Neumann  or  Robin.  The  chosen  solution  need  not  have  a  physical  meaning  since  verification 
(of  codes  or  of  calculations)  is  a  purely  mathematical  exercise.  But  choosing  a  physically 
realistic  manufactured  problem  which  has  a  closed  form  solution  offers  a  useful  advantage: 
It  exercises  each  term  involved  in  the  PDEs  in  a  manner  similar  to  that  of  a  real  problem 
so  that  similar  difficulties  in  the  solution  and  error  estimation  processes  will  arise. 

Once  the  manufactured  solution  has  been  constructed  and  the  source  terms  determined 
for  the  set  of  equations  to  be  verified,  code  verification  can  take  place  on  any  grid  in  the 
domain  covered  by  the  MMS.  By  verifying  the  code  on  increasing  grid  resolutions,  we  may 
show  that  as  grid  resolution  is  increased,  the  solution  becomes  more  accurate. 

2.4.2  Validation 

Once  a  CFD  code  has  been  verified,  it  is  important  to  validate  the  code’s  application 
to  certain  problems.  A  major  goal  of  this  thesis  is  to  validate  the  Spalart-Allmaras  model 
for  strand  grids.  Validation  of  this  model  within  the  strand  grids  framework  is  important 
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as  many  new  aspects  of  simulating  turbulence  are  investigated,  primarily  different  strand 
grid  meshing  strategies  to  best  resolve  the  flow  field  around  the  trailing  edge  of  an  airfoil 
and  to  capture  the  wake. 

Validation  may  be  described  as  “solving  the  right  equations.”  It  is  not  possible  to  vali¬ 
date  the  entire  CFD  code.  One  can  only  validate  the  code  for  a  specific  range  of  applications 
for  which  there  is  experimental  data.  Thus  one  validates  a  model  or  simulation.  Applying 
the  code  to  flows  beyond  the  region  of  validity  is  thus  termed  prediction.  A  more  precise 
definition  of  validation  in  this  scenario  is  given  as: 

The  process  of  determining  the  degree  to  which  a  model  is  an  accurate  representation  of  the 
real  world  from  the  perspective  of  the  intended  uses  of  the  model.  [21] 

Validation  examines  if  the  conceptual  models,  computational  models  as  implemented 
into  the  CFD  code,  and  computational  simulation  agree  with  real  world  observations.  The 
strategy  is  to  identify  and  quantify  error  and  uncertainty  through  comparison  of  simulation 
results  with  experimental  data.  The  experiment  data  sets  themselves  will  contain  bias  errors 
and  random  errors  which  must  be  properly  quantified  and  documented  as  part  of  the  data 
set.  The  accuracy  required  in  the  validation  activities  is  dependent  on  the  application,  and 
so,  the  validation  should  be  flexible  to  allow  various  levels  of  accuracy. 
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Chapter  3 

Governing  Equations 

This  chapter  describes  the  necessary  equations  needed  to  model  turbulent  flow  in  the 
strand  grids  framework.  The  general  equations  of  motion  for  a  fluid  will  be  described  first. 
How  these  equations  may  be  altered  to  accommodate  turbulence  modeling  will  follow.  The 
Spalart-Allmaras  model  will  then  be  described,  with  details  of  how  these  were  altered  for 
CFD  purposes. 


3.1  Fluid  Conservation  Laws 

To  accurately  model  laminar  two  dimensional  flow  in  CFD,  four  equations  must  be 
solved;  mass,  x-momentum,  y-momentum  and  energy.  These  equations  in  index  notation 
are  shown  below  in  Strong  Conservation  Law  (SCL)  form  in  equations  3.1  to  3.3. 


dp  d(puj) 
dt  '  dxj 


=  o 


d{pUj)  iHinijUj  •  !  d((Tjj)  .  _  q 

dt  '  dxj  '  dxj  '  r 


d(pe)  ,  djphuj ) 
dt  '  dxj 


_  dqj_ 

dxj  1  dxj 


pgjuj  =  0 


(3.1) 

(3.2) 

(3.3) 


Here  we  shall  define  p  as  pressure,  (ji  as  the  itli  component  of  the  body  force,  e  as  the  total 
energy  (internal  and  kinetic)  per  unit  mass  and  h  as  enthalpy;  which  may  also  be  defined 
as: 


h  =  e  +  2  (3.4) 

The  variable  qj  is  the  j-th  component  of  the  heat  flux  vector.  We  can  relate  the  heat  flux 
vector  to  the  temperature  gradient  through  Fourier’s  law  of  heat  conduction: 
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dT 

(l3  =  ~Kqlk 7 

(3.5) 

K  _  Cpfi  _  t Bjx 
'N  ~  Pr  -  (7_i  )pr 

R  is  the  gas  constant  for  the  fluid  in  use,  7  is  the  ratio  of  specific  heats,  Pr  is  the  Prandtl 
number  and  nq  is  the  coefficient  of  thermal  conductivity. 

Up  to  this  point,  the  only  assumptions  we  have  made  are  that  the  fluid  is  a  continuum, 
in  which  we  neglect  the  true  molecular  nature  of  matter.  We  will  make  two  other  assump¬ 
tions  about  the  equations  of  motion.  First,  we  will  neglect  body  forces,  such  as  gravity. 
This  is  a  good  approximation  for  isothermal  fluids  with  low  density,  such  as  air.  Second, 
we  will  assume  that  the  fluid  is  “Newtonian,”  which  means  that  the  rate  of  strain  of  the 
fluid  is  linearly  proportional  to  the  stress.  The  classic  Newtonian  stress  tensor  (ay/)  has  the 
form  [35]: 


aij  ~ 


2  duk  \ 

3  dxk  ' 


(3.6) 


From  a  CFD  perspective,  these  equations  may  be  written  in  a  more  compact  form  as: 


dQ  dFL_dJ!  = 

dr  '  dxj  dxj 


The  vectors  Q,  Fj  and  F?  are  defined  as: 


(3.7) 


(A 

(  PU3  ^ 

(  «  1 

Q  = 

put 

\pe) 

.  Fi  = 

pu.iUj  +  Sijp 

\  rjhu3  ) 

TTV  _ 

’  3 

^7 

<7 

1 

(3.8) 


Here,  Q  is  the  vector  of  conserved  variables.  These  are  the  dependent  variables  that  are 
solved.  Fj  is  known  as  the  inviscid  flux  vector,  and  FJ  is  known  as  the  viscous  flux  vector. 
If  we  set  Fj  =  0,  we  obtain  the  Euler  equations,  which  is  a  model  for  “inviscid  flow.”  The 
Euler  equations  provide  a  useful  “first  cut”  approximation. 
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Using  ordinary  CFD  techniques,  i.e.  not  DNS  or  LES,  the  solution  to  these  equations 
provide  an  accurate  representation  of  the  physical  problem  in  laminar  flow;  however,  if  the 
Reynolds  number  is  high  enough  such  that  the  flow  becomes  turbulent,  then  the  laminar 
solution  will  no  longer  be  accurate  or  appropriate.  Thus  to  accurately  model  turbulent  flow, 
these  equations  of  motions  must  be  modified  to  account  for  turbulence. 


3.1.1  Farve- Averaged  Navier-Stokes  Equations 

In  this  section  we  aim  to  look  at  both  the  Reynold-Averaged  Navier-Stokes  (RANS) 
equations  and  the  Farve- Averaged  Navier-Stokes  equations  (FANS).  Generally  the  RANS 
equations  are  appropriate  for  most  applications,  and  the  FANS  equations  used  only  when 
encountering  highly  compressible  or  hypersonic  flow.  As  the  strand-Cartesian  approach  is 
designed  for  aerodynamic  flows,  compressible  flows  will  be  encountered,  thus  we  use  FANS 
equations.  RANS  turbulence  models  are  equally  applicable  to  FANS  equations. 

Before  we  can  modify  the  equations  of  motion,  some  notation  must  first  be  introduced. 
We  introduce  classic  time  averaging  (Reynolds  averaging)  and  density  weighted  time  aver- 
gaing  (Farve  averaging).  Both  systems  are  broken  down  into  mean  and  fluctuations.  The 
mean  is  donated  by  either  an  overbar  or  tilde,  and  fluctuations  by  either  a  prime  or  double 
prime.  Reynolds  averaging  for  some  variable  <J>  is  defined  as: 


^  =  T  It 

$  =  $  +  <R 

Farve  averaging  for  some  variable  <F  is  defined  as: 


(3.9) 


| >  = 

p 

$  =  $  +  $" 


(3.10) 


It  is  important  to  note  that  the  mean  of  the  Reynolds  fluctuations  is  zero,  however,  this  is 
does  not  hold  true  for  Farve  fluctuations. 
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<!>'  =  0 
W7  /  0 


(3.11) 


Applying  both  Reynolds  and  Farve  averaging  to  the  equations  of  motion,  there  now 
exists  a  new  set  of  equations  which  account  for  turbulent  fluctuations.  Oliver  [36]  describes 
this  in  great  detail.  Taking  the  final  result  we  obatin  an  open  set  of  equations  that  is 
commonly  know  as  the  Farve- Averaged  Navier-Stokes  (FANS)  equations.  These  equations 
in  SCL  form  with  our  prior  assumptions  have  the  form  of: 


&p  I  9(pUj) 

dt  '  dxj 


=  0 


(3.12) 


8{pUi)  , 

dt  ‘+‘ 


d(pujui-pSij)  Q  / 


dxj 


+  Qx.XVii 


pu"u")  =  0 


(3.13) 


9(pE)  ,  d(pHuj) 
dt  '  dxj 


-  pu"u'-)}  -  {qj  +  pu”h")  +  £.{a, 


_d_t— 

dxj 


d 

dxj 


'V 


U ■  — 


2  PUjU'iu'i 


=  0 
(3.14) 


Immeadiately  we  can  identify  the  effects  of  turbulence  in  these  equations:  the  Reynolds 
stress  tensor  — pu'[u ”  and  the  Reynolds  heat  flux  pu'-h" .  Additionally,  we  see  two  terms  that 
refer  to  turbulent  transport  and  work:  transport  of  turbulent  kinetic  energy  by  turbulent 
velocity  fluctuations  —  ^ pu"u'(u "  and  work  done  by  the  viscous  stress  due  to  turbulent 
velocity  fluctuations  (?iju". 

We  define  some  turbulent  terms.  Total  turbulent  energy  is  defined  as  E  (total  energy 
and  turbulent  kinetic  energy),  H  is  the  total  turbulent  enthalpy,  k  is  the  turbulent  kinetic 
energy,  pt  is  the  turbulent  eddy  viscosity  and  qxj  is  the  turbulent  heat  flux.  We  also  define 
the  transport  and  work  of  turbulence  [36].  These  relationships  are  shown  below: 


E  =  e  +  k 


(3.15) 


H  =  h  +  k 


(3.16) 
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(3.17) 


d 

dxj 


{(Ji.ju'l 


(P  +  \ih) 


dk 

dxj 


(3.18) 


The  turbulent  Prandtl  number  is  Prx,  and  often  has  a  value  of  0.9  for  air.  For  comparison, 
the  Prandtl  number  for  air  often  has  a  value  of  0.75.  This  form  of  the  Navier-Stokes  can  be 
manipulated  into  a  vector  form,  resulting  in  [37]: 


(  \ 

(  \ 

p 

pUj 

Q  = 

put 

,  Fi  = 

piliUj  +  dijP 

\pE) 

\  PHui  ) 

pv  = 


Oij  -  pvffu'j 


\Ui(<Tij  -  pu'fu ")  +  (m  +  \pt)^ 


\ 


(Qj  +  QTj)J 

(3.19) 


Turbulence  models  are  based  on  solving  variables  so  that  the  Reynolds  stress  tensor 
may  be  given  a  value.  A  note  should  be  made,  that  from  this  point  on,  the  overbar  and 
tilde  on  singular  mean  variable  terms  shall  be  implied.  Fluctuating  terms  shall  still  retain 
their  prime  or  primes.  Many  turbulence  methods  utilize  the  Boussinesq  Approximation  in 
order  to  solve  for  this  term. 


3.1.2  Boussinesq  Approximation 

The  Boussinesq  Approximation  was  proposed  by  Boussinesq  in  1877  [38]  and  hypoth¬ 
esized  that  Reynolds  stresses  might  be  proportional  to  the  mean  rate  of  deformation.  This 
hypothesis  introduces  the  concept  of  eddy  viscosity  (//*).  Eddy  viscosity  describes  the  dis¬ 
sipation  of  momentum  caused  by  turbulent  eddies,  and  is  analogous  to  the  dissipation  of 
momentum  at  the  molecular  level  caused  by  a  fluid’s  viscosity. 


Tij  —  pu'fu'j  -  +  qx\  3dxkk^j)  \pk5, 


2  dub 


(3.20) 


The  turbulent  kinetic  energy  is  symbolized  by  k,  which  is  a  base  concept  for  many 
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turbulence  models,  most  notably  two-equation  models.  The  disadvantage  of  the  Boussinesq 
hypothesis  as  presented  is  that  it  assumes  pt  is  an  isotropic  scalar  quantity,  which  is  not 
strictly  true.  The  Reynolds  stress  tensor  shall  be  notated  by  from  this  point  on  to  keep 
notation  compact.  We  reach  an  equation  with  the  final  form  of: 


(  \ 

(  \ 

p 

puj 

Q  = 

put 

.  FJ  = 

pupij  +  SijP 

\PE  ) 

\  PHuJ  / 

rpv 

3 


&ij  T 


\ui(<Tij  +  Tij )  +  0  +  \pt)^ 


\ 


(Qj  +VTj)J 

(3.21) 


In  order  to  reduce  notation,  from  this  point  we  shall  let  +  rtj  and  qj  =  q3  +  qxj. 


3.2  Spalart-Allmaras  Turbulence  Closure 

In  Chapter  2,  the  general  form  of  the  Spalart-Allmaras  equation  was  shown  and  dis¬ 
cussed.  In  this  section,  more  precise  definitions  are  given  for  the  SA  model.  The  original 
Spalart-Allmaras  equation  has  the  form  of  [19]: 


pW,  =  pC„i(  1  -  MSP  - 1 /*Ww  -  2](|)2  +  ^llj((p  +  pp)£-)  +  cw>HI|] 

(3.22) 

Since  the  arrival  of  the  original  Spalart-Allmaras  model,  a  number  of  modifications  to 
the  model  have  been  made,  creating  new  and  alternate  SA  models.  Some  of  these  are  simple 
changes,  such  as  setting  trip  terms  to  zero.  Other  changes  are  more  drastic  such  as  those 
proposed  by  Catris  [39].  The  emergence  of  all  these  modified  SA  models  has  created  some 
confusion  when  selecting  which  form  of  the  SA  model  to  use. 

One  of  the  reasons  that  SA  models  were  altered  was  because  situations  arise  on  coarse 
grids  and  transient  states  where  the  turbulent  solution  may  become  negative.  This  is  often 
encountered  at  the  edge  of  the  boundary  layers  and  wakes,  where  the  turbulence  solution  is 
characterized  by  ramp  solutions  that  transition  to  constant  outer/freestream  regions  over 
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a  short  region.  The  rapid  transition  from  large  inner  to  relatively  small  outer  levels  can 
result  in  undershoots  for  discrete  solutions.  These  undershoots  may  cross  zero,  resulting 
in  a  negative  solution  values.  Previous  methods  of  eliminating  negative  solution  values  has 
been  to  clip  updates. 

Recently  Allmaras  has  added  some  much  needed  clarity  to  the  selection  of  SA  models, 
by  proposing  a  new  “negative  model”  [40],  which  is  to  be  used  in  conjunction  with  the 
original  SA  model  (or  another  SA  model  if  the  user  so  desires).  The  negative  SA  model 
corrects  these  discrete  undershoots  that  cross  zero,  and  gives  the  SA  model  robustness.  In 
this  work,  the  original  SA  model  (with  trip  terms)  shall  be  used  when  v  >  0.  When  z>  <  0  the 
negative  model  shall  be  utilized.  Allmaras’s  most  recent  work  makes  a  minor  modification 
to  the  original  proposed  SA  model.  This  modification  accommodates  compressible  flow 
by  not  holding  the  density  constant  in  the  turbulent  diffusion  coefficient  77  (see  below). 
By  not  holding  the  density  constant  the  SA  model  becomes  more  accurate  and  robust  in 
compressible  flow.  This  extra  term  is  absorbed  into  T  (see  below). 

To  make  any  SA  model  usable  in  CFD,  some  minor  alterations  must  be  made.  The 
first  modification  we  must  make  is  to  use  the  chain  rule  to  make  the  LHS  of  the  SA  equation 
comply  with  SCL  form.  This  modification  yields: 


Du  _  d(pu)  ,  d(puju) 
^  Dt  dt  '  dxj 


This  lone  modifications  yields  the  combined  standard-negative  model  below: 


(3.23) 


d(pu)  ,  d(puju)  _  g  (n  dv  \  I  p 
dt  '  dxj  dxj  '  <7  dxj '  ' 


D  +  T 


(3.24) 


The  source  terms  P,  D,  and  T  are  defined  explicitly  below.  P ,  the  rate  of  production 
of  v,  is  simply  known  as  the  production  term. 


v  >  0 


P=  < 


pCb  1(1  -  ft2)Sv, 

pCbl{  1  —  Ci3)flz>, 


V  <  0 


(3.25) 
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The  term  D  is  defined  as  the  rate  of  dissipation  of  D,  also  know  as  the  destruction  term. 


D  =  ( 


p[cuifu-^ft2](i)2,  ^>o 


-pCu,  l(|)2, 


v  <  0 


(3.26) 


It  is  important  to  note  that  when  v  <  0,  the  destruction  term  D  become  positive.  Below  T 
is  defined.  Commonly  known  as  the  transport  term,  it  is  a  term  for  the  transport  of  v  by 
turbulent  diffusion. 


rn  _  pCi2  dv  dv  _  {v+vfn)  dv  dp 

( j  dxk  dxk  a  dxk  dx ^ 


(3.27) 


The  turbulence  diffusion  coefficient  r/,  is  defined  as: 


r/  =  //  +  pvfn 


(3.28) 


fn 


1, 

Cnl+X3 
k  Cnl-X3  ’ 


V  >  0 

i>  <  0 


(3.29) 


Dynamic  eddy  viscosity  may  be  related  to  v  via  Eq  3.30.  When  u  is  negative,  the  eddy 
viscosity  is  set  to  zero,  and  u  itself  becomes  a  passive  scalar.  Using  this  relationship,  the 
Reynolds  stress  tensor  may  be  described  as  below: 


/R 


pvfv\,  v  >  0 

< 

0,  v  <  0 


(3.30) 


Tij  ~  Pvfvlidx]  +  dxi  3  sPkSij 


2  duk 


(3.31) 


The  turbulent  kinetic  energy  variable  is  not  calculated  in  the  Spalart-Allmaras  model,  and 
thus  the  last  term  is  often  ignored  when  estimating  Reynolds  stresses. 
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Wall  Damping  Functions  and  Closure  Coefficients 

Simple  eddy-viscosity  models  (EVM)  such  as  SA  rely  on  damping  functions  to  mimic 
the  attenuation  of  the  turbulence  near  the  wall.  These  damping  functions  depend  on  the 
local  turbulent  Re  number  or  wall  distance.  The  SA  model  contains  some  wall  damping 
functions,  such  as  the  one  mentioned  in  Eq  3.30.  The  wall  damping  function  fv i  has 
two  purposes.  At  the  wall,  the  damping  function  fv i  becomes  zero.  At  high  Reynolds 
numbers  the  function  tends  to  unity,  and  thus  kinematic  eddy  viscosity  v  simply  becomes 
the  kinematic  eddy  viscosity  z^. 


fv  1 


x 


3 


Xa+C«i 


X  = 


(3.32) 


The  wall  damping  function  fv 2  serves  as  similar  purpose  to  /,,  1 .  At  the  wall  the  damping 
function  becomes  unity,  and  at  high  Re  the  function  becomes  zero. 


fv  2  —  1  - 


\ 

i+xAi 


(3.33) 


Another  wall  function  used  is  the  non-dimensional  function  fw.  In  the  logarithmic  layer  of 
a  zero  pressure  gradient  boundary  layer  /„,  cs  J . 


fov 


l  +  g^.3 
96+Cl 


1 
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g  =  r  +  C^2(?’6  +  ?’) 


(3.34) 


r  —  u 

flK2d2 

d  is  the  distance  from  the  current  point  to  the  nearest  wall  and  tt  is  the  magnitude  of  the 
vorticity,  commonly  defined  by: 
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O  —  2  f  ^  j 


O . .  —  1  (  duj  _  duj 

&  —  o  I  a  .  i  i 
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ctei 


(3.35) 
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n  +  s, 

Q  Q(C2+C3S) 

(c3-2 c2)f2-S’ 


5  >  — C2^4 
5  <  — C2^4 


(3.36) 


O  ^  fv2 

‘  K2d2 ’ 


(3.37) 


It  is  important  to  note  that  12  must  not  be  allowed  to  reach  zero  in  order  to  avoid 
numerical  problems,  thus  this  particular  form  of  S  from  Oliver  [41]  is  used  to  avoid  this. 
The  closure  coefficients  suggested  by  Allmaras  [19,40]  are  as  follows: 


cb\  =  0.1355,  cb2  =  0.622,  a  =  §,  k  =  0.41 
Cuj 2  =  0.2,  3  =  2,  cv  1  =  7.1  Cti  =  1,  ct2  =  2,  ct3  =  1.2,  cti  =  0.5  (3.38) 


r  i  —  Cbl  -\- 


1+Cb2 

(T 


Final  SA  Form  Implemented  in  Strand  Grids 

Keeping  with  a  system  of  equation  notation,  the  Spalart- Allmaras  model  and  the  FANS 
equations  of  motion  may  be  written  as: 


dQ  ,  9El  _  91L  =  c 
dr  '  dxj  dxj  ‘~ 


(3.39) 


The  turbulent  kinetic  energy  k  has  been  ignored  in  the  FANS  energy  equation  (Eq.  3.14), 
since  SA  is  a  one-equation  model,  k  is  not  calculable.  The  vectors  Q,  Fj,  F?,  and  S  for  the 
standard-negative  SA  model  are  now  defined  as: 
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Chapter  4 

Numerical  Methods 

This  chapter  will  demonstrate  the  methods  in  which  the  governing  equations  provided 
in  Chapter  3  were  discretized  in  order  to  be  applied  to  strand  grids.  Note  that  in  the  strand- 
Cartesian  solver,  the  equations  presented  in  3.39  are  solved  as  coupled  equations.  While 
this  makes  the  implementation  and  solving  of  the  equations  more  expensive,  it  accounts  for 
all  cross-coupled  terms. 


4.1  Strand  Grids  Discretization  Method 

In  order  to  discretize  Eq.  3.39  using  the  finite  volume  method,  it  is  necessary  to  integrate 
over  the  volume  of  the  problem. 


iv 


8Q  dFj 
dt  dxj 


dFV 

-7T2-  =S)dV 
dxj 


(4.1) 


This  integral  can  be  simplified  using  the  divergence  theorem  (often  referred  to  as 
Gauss’s  theorem)  [42].  This  theorem  states  that  the  flux  through  the  boundary  is  equal  to 
the  sources  or  sinks  of  flux  in  the  closed  domain.  This  yields: 


iv 


dQ 

dt 


dV  + 


FjrijdA  — 


Fj’njdA  = 


SdV 


i v 


(4.2) 


It  is  important  to  maintain  that  integrals  with  respect  to  A  are  taken  at  the  boundaries 
of  the  control  volume.  Leibniz’s  theorem  states  that  for  a  fixed  control  volume,  the  volume 
is  not  a  function  of  time  and  the  time  integral  can  be  rewritten  [14]. 


QdV  + 
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(4.3) 
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Fig.  4.1:  A  sample  2D  mesh  with  cells  cl  and  c2.  k  represents  the  face  between  the  cells. 


To  simplify  the  derivation,  viscous  terms  are  assumed  to  be  negligible  and  the  discussion 
is  limited  to  the  inviscid  Euler  equations.  This  yields  the  following  equation: 


4-  !  QdV  +  [  FjrijdA  =  0 
dt  Jv  JA 


(4.4) 


The  derivation  thus  far  has  made  no  discrete  assumptions  and  has  been  entirely  contin¬ 
uous.  To  discretize  the  above  equation,  the  domain  is  split  into  cell-centered  control  volumes 
and  Q  is  assumed  to  vary  linearly  in  each  cell.  Q  over  the  entire  domain  is  piecewise  linear, 
which  creates  a  second-order-accurate  scheme. 

Discretely  defining  Q  as  linear  has  direct  effects  on  the  fluxes  on  the  boundary.  The 
fluxes  can  be  exactly  computed  on  the  boundaries  using  only  one  quadrature  point  (mid¬ 
point  rule).  A  sample  of  this  configuration  is  found  in  Fig.  4.1.  cl  and  c2  are  located  at 
the  center  of  the  cells  and  k  represents  the  quadrature  point  on  the  boundary.  The  inviscid 
fluxes  can  be  written  in  the  form  below: 


FjUjdA  =  /  FjrijdA^  ~  =  R(Q) 


(4.5) 
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J-k  represents  the  numeric  flux  at  the  quadrature  points  at  the  cell  interface  and  R  represents 
the  “residual.”  T  is  the  area  weighted  flux  given  below: 

T  =  FjUjA  (4.6) 

We  now  return  to  Eq  4.3  to  apply  this  formulation  to  include  viscous  and  source  terms. 
Changing  notion  slightly,  we  now  reside  in  cell  o  with  weighted  area  fluxes  between  i  cells 
around  it  in  pseudo-time,  yielding: 

d-^y0  +  E  -  E  -^  =  (4-7) 

i  i 

where  V0  is  the  volume  of  the  current  cell  o.  We  simplify  further  in  manner  similar  to  the 
above  method,  but  now  including  the  viscous  and  source  terms: 

RoiQ)  =  E  ^  -  E  Foi  -  SoVo  (4.8) 

i  i 

thus  getting: 

Cl^V0  +  ^(Q)  =  0  (4.9) 

This  gives  a  starting  point  to  now  apply  a  time-stepping  scheme  to  the  equation. 
Explicit  schemes  calculate  the  state  of  a  system  at  a  later  time  from  the  state  of  the  system 
at  the  current  time,  while  implicit  schemes  find  a  solution  by  solving  an  equation  involving 
both  the  current  state  of  the  system  and  the  later  one.  Mathematically,  if  Y (f)  is  the  current 
system  state  and  Y{t,  +  At)  is  the  state  at  the  later  time,  where  At  is  some  small  time  step, 
then,  for  an  explicit  method: 


Y(t  +  At)  =  F(Y(t)) 


(4.10) 


while  for  an  implicit  method  one  solves  an  equation: 
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G(Y(t),Y(t  +  At))  =  0  (4-11) 

in  order  to  find  Y(t  +  At). 

Using  k  as  the  current  pseudo-time  step,  and  k  +  1  as  the  next  pseudo-time  step  and 
applying  an  implicit  scheme  to  Eq  4.9  returns: 

nk- 1-1  _  Qk 

Vo  °  Ar  °  +  Ro(Qk+1)  =  0  (4.12) 

where: 

a  r>k 

Ro(Qk+1 )  =  Ro(Qk)  +  -qq(Qo+1  ~  Qo )  (4-13) 

applying  Eq  4.13  to  Eq  4.12  and  letting  A Qk  =  Qk+1  —  Q „  yields: 

^I+—j(AQk)  =  -Ro(Qk)  (4.14) 

where  I  is  the  indentity  matrix. 

Eq  4.14  is  the  general  form  that  the  strand  solver  uses.  It  is  in  this  equation  that 
we  shall  refer  to  the  RHS  and  LHS  from  here  on  out.  From  observation,  it  is  clear  that 
this  form  requires  the  calculation  of  the  residual  at  k  time  step  (RHS),  is  common  to  both 
explicit  and  implicit.  This  form  also  requires  the  calculation  of  the  derivative  of  Rk  with 
respect  to  Q  (referred  to  as  LHS  generally),  which  is  unique  to  an  implicit  scheme.  The 
next  section  describes  the  methods  used  to  calculate  the  residual  R0(Q )  and  the  residual 
Jacobian  Qjk-. 

4.2  Discretization  of  FANS  and  SA  Turbulence  Model  System 

This  section  directly  follows  from  the  previous  and  describes  the  necessary  discretization 
methods  required  for  the  strand  solver  in  an  implicit  scheme.  We  shall  break  R0(Qk )  back 
into  its  original  form  of  and  S0V0 
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4.2.1  Inviscid  Flux 

Returning  to  the  notation  formed  in  Eq  4.4  and  4.5,  J7/,.  is  defined  by  the  following: 

A  =  l  C HQl )  +  HQr))  -  Dk(Qn,  Ql)  (4.15) 

where  Dl(Qr,Ql)  is  the  artificial  dissipation,  Ql  and  Qr  represent  the  reconstructed 
terms. 

The  stability  of  Eq.  4.15  is  dependent  on  the  choice  of  the  artificial  dissipation  term. 
Without  dissipation,  the  numerical  flux  becomes  immediately  unstable.  Adding  a  proper 
dissipation  model  makes  Eq.  4.15  local  extremum  diminishing  (LED).  LED  schemes  require 
that  local  maxima  decrease  and  local  minima  increase.  LED  schemes  ensure  that  the  local 
extrema  are  bounded  and  prevent  them  from  diverging. 

Our  attention  now  turns  to  calculating  the  inviscid  part  of  To  calculate  the  LHS 
of  inviscid  FANS-SA  system  of  equations,  we  let: 


14 


dZ 

dQ 


i?|A|  R~l 


where  A  are  the  known  eigenvalues  defined  as: 


(4.16) 


A  =  diag(u',  u',  v!  +  c,  v!  —  c,  u ')  (4-17) 

We  define  c  as  the  speed  of  sound,  and  v!  =  un  —  vn.  We  now  have  un  defined  as: 


Axu-\-Avv 

Un  —  p|  —  W'xU'  HyV 

Knowing  A  and  A,  we  can  solve  for  R  and  i?-1. 


(4.18) 


4.2.2  Viscous  Flux 

The  viscous  portion  of  R %  and  must  now  must  be  calculated.  We  define  the 
weighted  area  viscous  flux  Fv  as: 
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F!  =  A-  Fv  =  AXFV  +  AyGv  =  Bvxf %  +  B'yf  (4. 19) 

Bxj  and  Byj  are  matrices  containing  the  constants  that  lie  in  front  of  their  respective 
differentiated  primitive  variables.  This  form  now  provides  a  format  to  solve  the  viscous 
portion  of  the  residual. 

0  0  0  0  0 

0  |  Ax/j  Ay/j,  0  0 

Bxf  =  0  -§  Ayn  Axn  0  0  |  (4.20) 

0  2/i(Axii—^)  n(Axv  +  Ayu)  Axk  0 

0  0  0  0  A ^ 

0  0  0  0  0 

0  Ay/j,  —  \AX  0  0 

B'yf  =  0  Axti  \Ayix  0  0  |  (4.21) 

0  n(Ayu  +  Axv)  2fi(Ayv  —  ^r)  Ayk  0 
0  0  0  0  Ayl 

The  viscous  flux  Jacobian  needs  to  now  be  addressed  to  satisfy  the  residual  Jacobian. 
We  wish  to  calculate  the  differential  of  the  viscous  flux  with  respect  to  Q  and  also  with 
respect  to  the  gradient  of  Q. 


at  a  fac e(/)  :  F>  =  Fv(Qf,  Qpf,  Qpf)  (4.22) 

Thus  we  have: 


8TV  _  8TV  9Qf  .  dTv  9®xf  dQp  ,  8TV  dQyf  dQp 
8Q  dQf  8Q  '  8Qpf  dQp  8Q  +  8Qpf  8QP  8Q 


(4.23) 


We  can  simplify  this  greatly  if  we  let  =  rp,,  =  Bxj,  =  Byj  and  =  Cj. 
Recall  that  Bxj  and  Byj  were  defined  previously. 
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a tv 

dQ 


S~1V  dQ  f 

°/  dQ 


+  [B 


v  a_Sk 

Xf  dQp 


+  B 


dQpf 

V  _^yf 

yf  dQp 


]rp 


(4.24) 


The  face  reconstruction  terms  must  then  be  linearized.  If  we  go  from  node  a  to  b 
through  a  face  with  nodes  1  and  2  on  each  corner  we  may  define  the  linearization  of  the 
reconstruction  terms  on  fine  levels  as: 


Cn  C 12 

C*21  C*22 


axaqp 


12 


1  V 


AyAQ 


(4.25) 


ab  , 


Choose  Qf  =  |( Qa  +  Qb),  thus  we  may  generalize  the  dehnition  as: 


DXAQP 


ab 


\yDyAQ 


(4.26) 


ab  > 


where  Dx  and  Dy  are  the  constants  that  change  with  respect  to  multi-grid  level  (either  fine 
or  course).  Putting  these  changes  together  produces: 


where: 


dZl 

dQ 


Cvf 


{h  +  ih)  +  B°f{  rf-rf) 


(4.27) 


Bf  =  DxBlf  +  DyBls 

A  final  form  of  the  viscous  Jacobian  can  be  written  as: 


dFv 

dQ 


{C}+BvfTp)b  +  {%~BvfTp)a 


(4.28) 


(4.29) 


4.2.3  Source  Terms 

Calculating  the  residual  portion  of  the  source  term  is  easy,  as  it  is  simply  the  source 
term  at  the  cell  multiplied  by  the  volume  of  the  cell.  Calculating  the  residual  Jacobian 
portion  of  the  source  term  however  is  a  complicated  process.  Similar  to  the  viscous  flux 
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Jacobian,  we  are  required  to  calculate  the  differential  of  the  source  term  with  respect  to  Q 
and  the  gradient  of  Q. 


S  =  p-D  +  T  =  P(Q,  VQ)  -  D(Q,  VQ)  +  T(Q,  VQ)  (4.30) 


An  example  of  how  the  source  terms  are  linearized  is  shown  below.  Source  terms  below 
are  differentiated  with  respect  to  Q;  P(Q),  D(Q ),  and  T(Q).  Note  that  all  these  derivatives 
have  been  greatly  simplified,  and  more  in-depth  derivatives  have  not  been  presented. 


dP 

dqi 


(  Cbi(§!j:pi>+^S)(l-ft2)+pi>SCbl^,  v>0 
cbl(i-Ct3)n^,  d  <  o 


(4.31) 


dD 

dqi 


d2  '  dqi 

Cwl  (  dpv2  - 
d2  1  dqi  ' 


Cui  ( dpv2  f  ,  dfti  ~2  \ 


+  Hf^2) -  +  itp"2),  s>o 


K,2d2  '  dqi 


dqi 


v  <  0 


(4.32) 


dT  Cb2  dv  dv  dp  _  1  dv  dp  r  dv  £  _i_  dfn  ~  ,  dv  \ 

dqi  ~  cr  dxk  dxk  dqi  a  dxk  dxk  1  dqi  ■’n  dqi  'dqi' 


(4.33) 


Below  we  show  the  source  terms  with  derivatives  with  respect  to  the  gradient  of  Q;  P(VQ), 
P>(VQ),  and  P(VQ). 
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dqXi 


pCbi(l  -  /t2)£^r,  v>0 

pCbl(l  -  Ct3)i> r<0 


dD 

dqXi 


I  pCuil 


(v\2  dfcj 
V  d  )  dqxi  ’ 


v  >  0 


0, 


v  <  0 


(4.34) 


(4.35) 


dT  _  pCb 2  d  i  dv  dv  \  _  (v+vfn)  d  r  dv  dp  \ 
dqxi  —  <y  dqxi  V  dxk  dxk  >  a  dqxi  V  dxk  dxk  ' 


(4.36) 


4.3  Boundary  Conditions  and  Boundary  Discretization 

Boundaries  for  strand  grids  are  very  similar  to  the  internal  discretization;  however, 
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unlike  internal  discretization  the  flux  on  the  wall  is  now  the  actual  problem  boundary 
instead  of  a  neighboring  cell.  Strands  use  only  a  cell-centered  structure.  Since  known 
values  in  cell-centered  paradigms  are  only  known  at  the  centroid  of  the  cell,  implementing 
equations  at  the  boundaries  is  difficult.  To  apply  equations  at  the  boundary,  a  “ghost-node” 
is  placed  at  the  quadrature  point  of  the  boundary  cell  (6).  The  values  from  the  interior 
of  the  domain  are  extrapolated  in  a  similar  fashion  to  the  interior  reconstruction  to  the 
boundary  (e).  Fig.  4.2  shows  a  simple  one  dimensional  example  of  this. 


Qe 


Qi  4“  (%e 


(4.37) 


Fig.  4.2:  A  sample  ID  cell  centered  boundary  with  extrapolated  point  and  ghost  node. 

The  general  form  of  Dirichlet  boundary  conditions  consists  of  B(QS )  =  b,  where  m 
number  of  equations  must  be  satisfied  for  each  condition.  Allmaras  [43]  and  Katz  [35] 
describe  this  in  detail  further.  The  boundary  fluxes  can  be  presented  in  the  form  of  a 
boundary  residual,  =  B(QS )  —  b,  similar  to  interior  fluxes.  Here  b  is  the  constant  values 
found  on  the  boundary,  and  B(QS )  some  auxiliary  variables  used  for  convenience,  primarily 
the  primitive  variables. 

In  turbulent  flow,  the  SA  model  requires  [19,40]  that  the  farfield  and  freestream  bound¬ 


ary  conditions  be  defined  as: 
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V  far  field  —  (3  t  5)^00  (4.38) 

and  the  wall  boundary  conditions  be  defined  as: 


b'wall  =  0  (4.39) 

Combining  the  above  SA  boundary  conditions  with  boundary  conditions  given  by  Katz 
[35],  we  may  obtain  the  equations  for  the  boundary  residual  for  key  boundary  conditions; 
inflow,  outflow,  viscous  wall,  inviscid  wall,  and  pure  dirichlet  where  the  entire  state  is 
specified  and  fixed.  The  spec  values  refer  to  specified  user  values.  An  example  of  an  outflow 
that  specifies  a  pressure  has  a  boundary  residual  of: 


Rb 


Pb  Pspec 
U  —  Ue 

<  V  —  Ve  > 

Tb-Te 


=  0 


An  Inflow  that  specifies 


a  pressure  has  a  boundary  residual  of: 


(4.40) 


(  Pb  Pspec  1 


Ut 


Rb={ 


Un 


=  0 


Tb-Te 


(4.41) 


^  b'gpec  J 


Here,  ub  =  u  •  t  is  the  velocity  tangential  to  the  inflow  plane,  along  the  tangential  vector, 
t  =  (—ny,nx).  In  a  similar  fashion,  un  =  u  •  n  is  the  velocity  normal  to  the  inflow  plane 
along  the  normal  vector,  n  =  ( nx,ny ).  The  viscous  wall  boundary  residual  is  given  as: 
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Pb  ~Pe 
Ub  +  Ue 
Rb  =  {  Vb  +  Ve 

Tb-Te 

^b  +  b'e 

An  inviscid  wall  boundary  residual  is  given  as: 


>  =0 


(4.42) 


P~  Pe 
Ut 


Rb=  { 


V* n 


>  =o 


(4.43) 


pe  -  pee 


[pp-  ppe  ) 


The  inviscid  wall  auxiliary  variables  Qs  use  the  standard  conserved  variables  seen  in  the 
equations  of  motion  instead  of  the  primitive  variables.  The  Dirichlet  boundary  residual  is 
given  as: 


Rb 


P  Pspec 
U  Uspec 
*  V  VSpec  > 

T  —  T 

J  ■L  spec 
b'  l^spec 


=  0 


(4.44) 


4.4  Method  of  Manufactured  Solutions 

To  start  the  Method  of  Manufactured  Solutions,  we  must  first  give  values  to  the  con¬ 
served  variables.  The  values  generally  consist  of  some  positive  smooth  trigonometric  func¬ 
tion.  These  terms  for  the  conserved  variables  shall  be  called  the  source  variables  to  avoid 
confusion.  An  example  of  how  this  is  done  for  density  is  shown  below: 
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p  =  B(  1  +  axsin(bx  -  cxx)  +  aysin(by  -  cyy)  +  axysin(bxy  -  cxyxy ))  (4.45) 

Here  x  and  y  represent  Cartesian  coordinates  along  the  grid.  B  is  some  base  value 
defined  by  the  user  to  keep  the  values  in  the  solution  realistic.  ax,  ay,  axy,  bx .  by,  bxy,  cx, 
cy ,  and  cxy  all  represent  some  arbitrary  constants  chosen  to  keep  the  function  smooth  and 
positive.  This  type  of  function  is  done  for  all  conserved  variables. 

From  3.40  we  show  the  equation  used  for  density  in  the  fully  coupled  FANS  and  SA 
equation: 

f  +  =  0  (4.46) 

applying  MMS  to  this  equation  in  steady  state  (§j  =  0),  we  let  this  equation  equal  to  some 
known  source  term,  denoted  by  Sp. 


9(pUj)  _  n 

dxj  ~  JP 


(4.47) 


The  source  term  is  calculated  directly  by  differentiating  the  source  variables  and  ap¬ 
plying  the  constants  defined  above.  In  the  2D  system  for  p  we  declared  above,  we  have: 


Sp  —  Fx  +  Gy 

f.=P5  +  u%  <448> 

O,  =P%+v% 

So  in  this  example,  we  shall  differentiate  our  source  variables  u  and  p  from  Eq  4.45 
with  respect  to  x  and  y  to  find  values  for  Fx,  Gy,  and  ultimately  Sp.  With  a  value  for  Sp .  we 
now  feed  the  source  variables  into  the  strand  solver.  If  the  algorithms  have  been  properly 
implemented,  then: 


d(puj) 

dxj 


Sp  =  error 


(4.49) 


51 


38 

As  the  strand  solver  is  not  exact  and  makes  some  assumptions  and  truncations,  there 
will  be  some  small  difference  {error)  between  the  source  term  and  the  solution.  We  perform 
this  test  for  increasing  grid  resolutions  to  show  that  as  the  number  of  cells  increase  for  a 
fixed  mesh  size,  the  accuracy  of  the  solution  improves  by  an  order  of  two.  The  results  of 
this  test  is  shown  in  Chapter  5. 
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Chapter  5 

Results  and  Discussion 

This  section  will  present  and  discuss  the  results  from  the  verification  and  validation  of 
the  Spalart-Allmaras  turbulence  model. 

5.1  Verification 

First,  we  perform  grid  refinement  studies  using  the  method  of  manufactured  solutions 
[22].  This  study  ensures  to  the  extent  possible  that  we  have  not  made  coding  mistakes  in 
the  strand  grid  discretization. 

Five  levels  of  grid  refinement  are  used  for  the  study,  with  32,  64,  128,  256,  and  512 
cells  along  each  side  of  the  square,  respectively.  The  manufactured  solution  is  chosen  based 
on  smooth  trigonometric  functions  similar  to  previous  work  [44] .  An  example  of  how  these 
functions  were  set  up  and  applied  to  the  strand  solver  was  shown  in  Chapter  4. 


(a)  Density  manufactured  solution  used  in  grid  re-  (b)  RMS  density  error  for  manufactured  solution 
finement  study.  grid  refinement  study. 


Fig.  5.1:  MMS  grid  refinement  study  results. 
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The  manufactured  solution  for  density  is  shown  in  Figure  5.1(a).  Dirichlet  conditions 
enforcing  the  exact  solution  are  used  at  all  boundaries  for  this  study.  The  results  of  the  grid 
refinement  study  for  isotropic  strand  distributions  are  shown  in  Figure  5.1(b).  h  is  defined 
as  the  characteristic  cell  size;  the  strand  length  (1),  over  the  number  of  cells  per  strand. 
Second  order  convergence  was  achieved  for  all  conserved  variables. 

Work  et  al.  [12]  performed  grid  refinement  studies  using  MMS  on  multi-strand  and 
strand  smoothing  for  laminar  flow,  allowing  for  any  sensitivities  of  the  flow  solver  to  the 
different  meshing  strategies  to  be  identified.  MMS  was  set  up  in  the  same  manner  as  it 
was  in  this  work.  Three  sharp  corner  treatments  were  investigated:  no  smoothing  or  multi¬ 
strands  (Figure  5.2(a)),  extremely  smoothed  strands  (Figure  5.2(b)),  and  multi-strands 
(Figure  5.2(c)).  We  perform  the  refinement  study  for  both  isotropic  (Figure  5.3(a))  strand 
distributions,  as  well  as  anisotropic  (Figure  5.3(b))  strand  distributions  for  which  the  wall 
spacing  is  chosen  to  produce  an  aspect  ratio  of  103  for  the  first  cell  off  the  wall. 


(a)  No  smoothing. 


(b)  Extreme  smoothing.  (c)  Multi-strand. 


Fig.  5.2:  Various  meshing  strategies  used  for  manufactured  solution  grid  refinement  studies. 


The  results  of  the  grid  refinement  study  for  isotropic  and  anisotropic  strand  distribu¬ 
tions  are  shown  in  Figure  5.3.  For  both  isotropic  and  anisotropic  distributions,  no  special 
corner  treatment  (labeled  “No  Smoothing”)  results  in  loss  of  error  convergence.  This  is 
expected  because  without  smoothing,  the  grid  cells  are  not  reduced  in  size  properly  as  the 
grid  is  refined.  Clearly,  special  treatment  is  required  at  sharp  corners  to  achieve  the  design 
order  of  accuracy. 


RMS  Density  Error 
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(a)  Aspect  ratio  of  1  -  grid. 


(b)  Aspect  ratio  of  103  -  grid. 


1/h 

(d)  Aspect  ratio  of  103  -  error. 


Fig.  5.3:  RMS  density  error  for  manufactured  solution  grid  refinement  studies  using  various 
meshing  strategies  on  isotropic  and  stretched  strand  distributions. 
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Second  order  accuracy  is  obtained  for  both  smoothing  and  multi-strand  approaches, 
for  both  isotropic  and  anisotropic  strand  distributions.  The  multi-strand  approach  gives 
slightly  lower  errors  with  isotropic  strand  distributions,  perhaps  due  to  reduced  skewing 
of  the  grid  cells.  But  for  anisotropic  high  aspect  ratio  strand  distributions,  which  would 
be  used  for  turbulent  RANS  or  FANS  calculations,  smoothed  strands  produce  lower  error 
than  the  multi-strand  grids  by  about  a  factor  of  five.  The  explanation  for  this  likely  lies 
in  the  interaction  of  multi-strands  and  stretched  node  distributions,  as  shown  in  the  red 
boxes  of  Figure  5.2(c).  Multi-strands  tend  to  produce  highly  discontinuous  mesh  spacing  at 
different  orientations,  both  at  the  strand  ends,  as  well  as  the  surface.  Attempts  were  made 
to  cluster  surface  points  near  the  corners  to  avoid  the  thin  high  aspect  ratio  cells  bordering 
the  tiny  wedge  cells  at  the  sharp  corner.  However,  in  doing  this,  the  mesh  discontinuity  at 
the  strand  ends  became  more  severe.  It  appears  that  eliminating  these  mesh  discontinuities 
could  prove  difficult  which  is  unfortunate  since  they  seem  to  be  responsible  for  higher  errors 
as  well  as  some  convergence  issues. 

5.2  Validation 

Extensive  use  of  the  NASA  Langley  turbulence  modeling  resource  has  been  made  for 
these  cases.  This  resource  is  comprised  of  previous  experimental  results,  theoretical  results, 
and  two  independent  compressible  CFD  codes;  FUN3D  and  CFL3D  [3]. 

5.2.1  Flat  Plate  with  Zero  Pressure  Gradient 

MMS  allows  for  the  solver  to  be  verified  on  a  mathematical  level,  however,  there  may 
still  be  problems  with  the  code  when  running  a  real  case.  Thus  additional  verification  is 
required.  The  turbulent  flat  plate  validation  case  may  be  considered  as  verification  through 
validation  due  to  its  simplicity. 

The  turbulent  flat  plate  case  was  run  at  M  =  0.2,  at  a  Reynolds  number  of  Re  =  5  x  106 
based  on  length  “1”  of  the  grid.  The  body  reference  length  is  2  units.  The  following  plot 
shows  the  layout  of  the  simple  flat  plate  grids  used  for  this  study,  along  with  the  boundary 
conditions.  An  outflow  BC  was  used  instead  of  a  farfield  Riemann  BC. 
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Flat  Plate  Boundary  Conditions, 

M=0.2,  ReL  =  5  million  (L=  1  >,  Tref  =  540  R 


0.5 


0  - 


farfield  Riemann  BC 


-  Pt/Pref=  1.02828, 

Tt/T,ef=  '•008> 

1  quantity  from  interior 


P/P^I.D, - ► 

other  quantities 
from  interior 


symmetry 


adiabatic  solid  wall 


start  of  plate  at  x=0 

J  i  i  i  i  1  i  i  i  i  L. 


-0.5 


0.5 


1.5 


Fig.  5.4:  Flat  Plate  Dimensions  and  Boundary  Conditions  [3] 


Certain  inputs  and  slight  deviations  were  required  to  each  model  in  order  to  comply 
with  the  NASA  Langley  turbulence  resource  validation  case.  The  SA  model  only  required 

that  V far  field  —  3r/oc. 

The  grid  used  for  the  flat  plate  validation  is  the  137  x  97  grid  shown  in  Figure  5.5(a). 
The  plate  leading  edge  begins  at  x  =  0  and  extends  for  a  length  of  two.  A  short  inviscid 
wall  entry  way  beginning  at  x  =  —0.33  is  provided  to  allow  for  proper  inflow  conditions. 
Stagnation  temperature  and  pressure  are  specified  at  the  inflow,  and  static  pressure  is 
specified  at  the  outflow.  The  turbulent  viscosity  held  for  this  case  is  shown  in  Figure  5.5(b) 
which  has  been  scaled  by  a  factor  of  40  vertically  to  facilitate  visualization.  Streamwise 
velocity  and  turbulent  viscosity  profiles  are  shown  for  two  locations  downstream  on  the 
plate,  and  are  over  plotted  with  FUN3D  and  CFL3D  results  in  Figure  5.6. 

Note  that  good  agreement  is  obtained,  even  for  the  137  x  97  grid,  which  is  16  times 
more  coarse  than  the  FUN3D  and  CFL3D  results  shown  in  the  figure.  The  computed  drag 
coefficient,  which  is  entirely  due  to  skin  friction  for  this  case,  is  shown  in  Table  5.1,  along 
with  FUN3D  and  CFL3D  results  for  the  same  grid.  The  drag  coefficient  falls  within  the 
range  predicted  by  the  established  codes. 
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(a)  137x97  grid  (b)  Contours  of  nt/  po 

Fig.  5.5:  Grid  and  turbulent  viscosity  contours  for  flow  over  a  flat  plate  at  M  =  0.2  and 
Re  =  5  x  106. 

The  pv  field  for  this  case  is  shown  in  Figure  5.7  which  has  been  scaled  by  a  factor  of  40 
vertically  to  facilitate  visualization.  It  is  clear  that  v  reaches  a  negative  value  in  the  flow 
field.  This  plot  confirms  the  need  for  the  “negative”  model. 

5.2.2  Bump-in-channel 

The  bump-in-channel  case  is  different  from  the  simpler  flat  plate  verification  case  be¬ 
cause  it  involves  wall  curvature  and,  as  a  result,  pressure  gradients.  It  was  run  at  a  Mach 
number  of  M  =  0.2,  at  a  Reynolds  number  of  Re  =  3  x  106  based  on  length  “1”  of  the  grid. 
The  body  reference  length  is  1.5  units.  The  lower  wall  is  a  viscous-wall  bump  extending 
from  x  =  0  to  1.5  The  maximum  bump  height  is  Yq  =  0.05.  The  exact  definition  of  the 
bump  is: 


Table  5.1:  Comparison  of  computed  drag  coefficients  for  flow  over  a  flat  plate  at  M  =  0.2 
and  Re  =  5  x  106. 


cd 

Strand 

2.82287E-3 

FUN3D  (quads) 

2.84005E-3 

FUN3D  (triangles) 

2.80289E-3 

CFL3D 

2.86621E-3 
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(a)  x-velocity  profiles  at  x  =  0.97008. 
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(b)  x-velocity  profiles  at  x  =  1.90334. 
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Fig.  5.6:  Comparison  of  streamwise  velocity  and  turbulent  viscosity  profiles  for  flow  over  a 
flat  plate  at  M  =  0.2  and  Re  =  5  x  106. 


y=  < 


0.05(sm(vr^  -  (f  )))4,  0.3  <  x  <  1.2 


(5.1) 


0, 


0  <  x  <  0.3  and  1.2  <  x  <  1.5 


The  upstream  and  downstream  farfield  extends  25  units  from  the  viscous-wall,  with 
inviscid  plane  boundary  conditions  imposed  on  the  lower  wall  between  the  farfield  and  the 
solid  wall.  The  upper  boundary  is  a  distance  of  y  =  5.0  high.  It  is  taken  to  be  an  inviscid 
plane.  Figure  5.8  shows  the  layout  of  this  case,  along  with  the  boundary  conditions. 

The  grid  used  for  the  bump-in-channel  validation  is  the  353  x  161  grid  shown  in  Figure 
5.9(a)  and  a  close  up  is  shown  in  Figure  5.9(b).  Figures  5.9  and  5.9(c)  are  scaled  by 
a  factor  of  100  vertically  and  a  factor  of  16.667  horizontally  to  facilitate  visulization.  To 
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Fig.  5.7:  Contours  of  pD 
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Bump-in-channel  Boundary  Conditions, 
M=0.2,  ReL  =  3  million  (L=1 ),  Trel  =  540  R 
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(a)  Dimensions  and  Boundary  Conditions 


(b)  Bump  close-up 


Fig.  5.8:  Bump-in-channel  Dimensions  and  Boundary  Conditions  [3] 


avoid  the  internal  corner  issues  over  the  bump,  strands  were  set  vertically  only,  and  were  not 
smoothed  or  skewed.  Stagnation  temperature  and  pressure  are  specified  at  the  inflow,  and 
static  pressure  is  specified  at  the  outflow,  however  similar  input  conditions  will  yield  similar 
results.  The  turbulent  viscosity  field  for  this  case  is  shown  in  Figure  5.9(c).  Streamwise 
velocity  and  turbulent  viscosity  profiles  are  shown  for  two  locations  downstream  on  the 
bump,  and  are  over  plotted  with  FUN3D  and  CFL3D  results  in  Figure  5.10. 

Note  that  good  agreement  is  obtained  for  all  profiles,  even  for  the  353  x  161  grid,  which 
is  16  times  more  coarse  than  the  FUN3D  and  CFL3D  results  shown  in  the  figure.  Below  a 
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(a)  353x161  grid 


Fig.  5.9:  Grid  and  turbulent  viscosity  contours  for  flow  through  a  bump- in- channel  at 
M  =  0.2  and  Re  =  3  x  106. 


plot  of  the  surface  coefficient  of  pressure  and  friction  along  the  bump  is  shown,  and  is  over 
plotted  with  CFL3D  and  FUN3D  results  in  Figure  5.11. 

The  NASA  Langley  turbulence  resource  [3]  makes  some  useful  observations  about  the 
surface  skin  friction.  In  this  bump  case  the  surface  skin  friction  is  singular  (tends  toward 
infinity)  at  the  leading  edge.  The  finer  the  grid,  the  more  nearly  singular  the  local  behavior 
on  a  finite  grid.  There  is  also  locally  anomalous  behavior  in  Cf  at  the  back  end  of  the  bump 
wall  (at  x  =  1.5),  as  is  often  seen  in  CFD  solutions  near  trailing  edges  [45].  Figure  5.11 
shows  the  local  anomalous  behavior  that  exists  near  the  leading  edge  ( x  =  0)  due  to  singular 
behavior  of  the  solution,  as  well  as  near  the  trailing  edge  (x  =  1.5)  due  to  numerical 
influences.  These  behaviors  differ  for  the  three  codes,  and  result  in  small  local  deviations 
that  can  be  seen  when  zoomed  into  the  two  locations.  All  three  codes  are  seen  to  yield 
nearly  identical  results  over  most  of  the  bump  wall. 

The  computed  drag  coefficient,  is  shown  in  Table  5.2,  along  with  FUN3D  and  CFL3D 
results  for  the  same  grid.  The  drag  coefficient  falls  within  the  range  predicted  by  the 
established  codes.  It  must  be  emphasized  again  that  the  C d  results  in  Table  5.2  obtained 
by  FUN3D  and  CFL3D  were  run  on  very  fine  grids;  while  the  strand  code  ran  on  a  grid 
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(b)  x-velocity  profiles  at  x  =  1.2. 


Ik/v  o 


(c)  turbulent  viscosity  profile 


Fig.  5.10:  Comparison  of  streamwise  velocity  and  turbulent  viscosity  profiles  for  flow 
through  a  bump-in-channel  at  M  =  0.2  and  Re  =  3  x  106. 


(b)  Surface  coefficient  of  friction 


Fig.  5.11:  Bump-in-channel  surface  coefficient  of  pressure  and  friction 
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16  times  coarser.  In  the  NASA  Langley  turbulence  resource,  CFL3D  and  FUN3D  obtained 
higher  values  of  Cd  on  coarse  grids.  As  the  number  of  cells  in  the  grid  increased,  both  codes 
tended  asymptotically  towards  the  numbers  given  in  Table  5.2. 

5.2.3  NACA  0012  Airfoil 

With  the  SA  model  verified  for  the  flat  plate  and  bump-in-channel  cases,  we  may 
validate  the  flow  over  a  NACA  0012  airfoil  at  M  =  0.15  and  Re  =  6  x  106  at  various  angles 
of  attack.  To  reduce  the  effect  from  the  farfield  boundaries,  the  Cartesian  grid  is  extended 
500  chords  in  all  directions.  Riemann  invariant  boundary  conditions  are  imposed  at  the 
farfield.  Portions  of  the  Cartesian  grid  system  are  shown  in  Figure  5.12.  A  baseline  grid 
consisting  of  10  levels  of  refinement  and  16,348  cells  is  used,  shown  in  Figure  5.12(a). 

Additionally,  a  wake-capturing  Cartesian  grid  consisting  of  15  levels  of  refinement  and 
304,348  cells  was  used,  shown  in  Figure  5.12(b).  The  wake  capturing  grid  is  made  quite 
wide  and  fine  to  handle  high  angles  of  attack.  The  grid  count  could  be  significantly  reduced 
by  using  a  high-order  method  and  automatic  grid  adaptation,  though  these  avenues  are  not 
pursued.  Nonetheless,  the  wake-capturing  Cartesian  grid  is  still  quite  efficient,  due  to  the 
efficiency  of  the  block-structured  Cartesian  grid  code. 

The  strand  grid  for  this  case  consists  of  a  NACA  0012  airfoil  containing  320  surface 
nodes  and  64  cells  along  each  strand,  yielding  20,480  total  strand  cells.  Both  multi-strands 
and  strand  smoothing  are  tested,  with  and  without  the  wake-refined  Cartesian  grid.  The 
four  mesh  combinations  are  shown  in  Figure  5.13.  The  wake-refined  Cartesian  grids  are 
wider  than  would  normally  be  needed  and  could  be  tailored  to  specific  cases  through  more 
sophisticated  Cartesian  mesh  adaptation  strategies. 


Table  5.2:  Comparison  of  computed  drag  coefficients  for  flow  through  a  bump-in-channel 
at  M  =  0.2  and  Re  =  3  x  106. 


cd 

Strand 

3.57559E-3 

FUN3D 

3.56106E-3 

CFL3D 

3.57238E-3 
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(b)  Wake  capturing  Cartesian  grid. 


Fig.  5.12:  Grid  system  for  flow  over  a  NACA  0012  airfoil. 


The  surface  pressure  coefficients  for  a  =  0°,  10°,  15°  are  shown  in  Figure  5.14,  and 
compared  with  the  experimental  data  of  Gregory  and  O’Reilly  [46].  The  Gregory  data  is 
actually  taken  at  Re  =  3  x  106,  not  Re  =  6  x  106,  but  little  change  in  pressure  and  lift  is 
observed  between  the  two  Reynolds  numbers.  The  figure  only  displays  the  results  for  multi¬ 
strand  treatment  without  the  wake-refined  Cartesian  grid  because  all  methods  matched  the 
pressure  data  reasonably  well. 

Greater  discrepancies  in  the  methods  are  uncovered  by  observing  computed  lift,  and 
drag  values.  Figures  5.15  and  5.16  show  lift  coefficient  versus  angle  of  attack  and  drag  co¬ 
efficient  versus  lift  coefficient,  respectively,  along  with  the  corresponding  experimental  data 
of  Ladson  [47] .  The  lift  data  is  matched  reasonably  well  by  most  mesh  conhgurations,  with 


Table  5.3:  Comparison  of  computed  lift  and  drag  coefficients  for  flow  over  a  NACA  0012 
airfoil  at  M  =  0.15,  a  =  15°  Re  =  6  x  106.  Note  that  the  FUN3D  and  CFL3D  results  use 
a  much  finer  grid. 


Ci 

cd 

Multi-strand,  no  Cart,  refine  (320  points) 

1.5587 

0.02332 

Multi-strand,  Cart,  refine  (320  points) 

1.5457 

0.02284 

Smoothed  strand,  no  Cart,  refine  (320  points) 

1.6127 

0.02417 

Smoothed  strand,  Cart,  refine  (320  points) 

1.5578 

0.02189 

FUN3D  (513  points) 

1.5547 

0.02159 

CFL3D  (513  points) 

1.5461 

0.02124 

(c)  Smoothed  strands,  no  wake  grid.  (d)  Smoothed  strands,  wake  grid. 

Fig.  5.13:  Various  meshing  strategies  for  the  trailing  edge  of  a  NACA  0012  airfoil. 
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(a)  a  =  0° 


(b)  a  =  10° 


Fig.  5.14:  Pressure  coefficient  for  flow  over  a  NACA  0012  airfoil  at  M  =  0.15  and  Re  = 
6  x  106  at  various  angles  of  attack.  Only  multi-strand  results  with  no  wake  Cartesian 
refinement  are  shown.  Other  meshes  were  indistinguishable. 


the  exception  of  the  smoothed  strand  case  with  no  Cartesian  refinement,  which  significantly 
over  predicts  lift  for  a  =  15°.  This  case  also  produces  the  most  error  in  the  drag  for  all 
angles  of  attack.  Importantly,  the  smoothed  strands  actually  match  the  data  better  than 
any  other  configuration  with  the  use  of  the  wake-refined  Cartesian  grid.  The  multi-strand 
cases  improve  less  in  the  presence  of  the  wake-refined  Cartesian  grid.  For  clarity,  the  drag 
results  for  a  =  15°,  shown  in  Figure  5.16,  are  summarized  in  Table  5.3,  along  with  FUN3D 
and  CFL3D  results  using  a  much  finer  grid  (513  surface  nodes  instead  of  320).  Again,  it  is 
clear  that  the  smoothed  strands  produce  the  worst  results  of  any  case  without  the  wake- 
refined  grid,  but  the  best  results  of  any  case  with  the  wake-refined  grid.  The  multi-strand 
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Fig.  5.15:  Ci  vs.  a  compared  to  experiment  using  various  meshing  strategies  for  flow  over 
a  NACA  0012  airfoil  at  M  =  0.15  and  Re  =  6  x  106. 

configurations  lie  in  the  middle  and  are  less  affected  by  the  wake-refined  Cartesian  grid. 

The  above  observations  are  explained  by  examining  the  ability  of  the  various  mesh 
configurations  to  accurately  resolve  the  thin  wake  leaving  the  airfoil  trailing  edge.  Resolution 
of  the  wake  is  critical  for  accurately  computing  lift,  and  especially  drag,  for  this  case.  Field 
plots  of  velocity  magnitude  show  the  ability  of  the  strand  method  to  accurately  capture  the 
wake  for  various  angles  of  attack  in  Figure  5.17.  The  figures  on  the  left  use  smoothed  strands 
without  the  wake-refined  Cartesian  grids,  while  the  figures  on  the  right  use  smoothed  strands 
with  the  wake-refined  grids.  Clearly,  the  refined  Cartesian  grids  capture  the  wake  much 
better,  which  explains  the  improved  drag  results  in  Table  5.3.  But  there  is  a  more  subtle, 
yet  critical  observation  beyond  this.  In  the  multi-strand  case,  the  increased  refinement  at 
the  trailing  edge  provided  by  the  multi-strands  provides  reasonable  resolution  of  the  wake, 
even  without  the  refined  Cartesian  grid.  In  the  smoothed  case,  however,  the  resolution  is 
much  worse  at  the  trailing  edge  than  is  the  case  with  multi-strands,  which  explains  why 


67 


54 


0.025 

0.020 

0.015 

0.010 

0.005 


□  Ladson 

•  multi-strand,  no  Cart,  refine 

. . 

A 

;  +  multi-strand,  Cart,  refine 

f 

a  smoothed  strand,  no  Cart,  refine 

□ 

*  smoothed  strand,  Cart,  refine 

□ 

□ 

▲  □ 

*  □  n 

:  □ 

□ 

□ 

□ 

00  05  L0  L5 


Fig.  5.16:  Ci  vs.  Cd  compared  to  experiment  using  various  meshing  strategies  for  flow  over 
a  NACA  0012  airfoil  at  M  =  0.15  and  Re  =  6  x  106. 


smoothed  strands  are  more  reliant  on  the  Cartesian  grids  to  obtain  accurate  results.  Also 
notable  are  the  lack  of  mesh  discontinuities  in  the  smoothed  grid  which  leads  to  more 
accurate  results  than  the  multi-strand  approach. 

In  summary,  the  optimal  mesh  configuration  for  this  case  appears  to  be  smoothed 
strands  with  wake-refined  Cartesian  grids.  Use  of  the  smoothed  strands  avoids  mesh  dis¬ 
continuities  present  with  multi-strands,  although  the  accuracy  of  smoothed  strands  is  highly 
reliant  on  the  availability  of  Cartesian  refinement  for  proper  resolution  of  the  wake  and  accu¬ 
rate  drag  computation.  The  use  of  a  high-order  Cartesian  method,  such  as  in  the  Helios  [48], 
is  expected  to  dramatically  improve  wake  resolution  even  further,  as  well  as  resolve  other 
important  features,  such  as  vortical  shedding. 

5.2.4  NACA  4412  Trailing  Edge  Separation 

With  the  simpler  NACA  0012  airfoil  case  validated,  a  more  complex  airfoil,  in  the  form 
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(e)  no  wake  grid,  a  =  15°  (f)  wake  grid,  a  =  15° 


Fig.  5.17:  Field  plots  of  velocity  magnitude  for  various  cr,  showing  the  effect  of  Cartesian 
refinement  in  the  wake  region.  All  cases  use  smoothed  strands. 
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of  the  NACA  4412  airfoil,  can  now  be  tested.  This  case  aims  to  study  the  flow  separation 
at  the  trailing  edge  of  the  airfoil  and  flow  characteristics  at  the  wake. 

For  purposes  of  the  validation,  the  definition  of  the  airfoil  shape  is  slightly  altered  so 
that  the  airfoil  closes  at  chord  x  =  1  with  a  sharp  trailing  edge. 

Experimental  test  data  is  provided  by  Coles  and  Wadcock  [49, 50]  Flowfield  character¬ 
istics  were  measured  with  a  flying  hot-wire  for  the  airfoil  at  13.87  degrees  angle  of  attack. 
A  Mach  number  of  M  =  0.09  at  a  Reynolds  number  of  Re  =  1.52  x  106  per  airfoil  chord  was 
used.  Both  the  upper  and  lower  boundary  layers  were  tripped  in  the  experiment  (2.5c  on 
the  upper  surface  and  10.3c  on  the  lower  surface).  Like  the  other  validation  cases  performed 
so  far,  FUN3D  and  CFL3D  will  be  used  for  comparative  purposes.  Note  that  FUN3D  and 
CFL3D  have  grids  with  a  farfield  outer  boundary  extending  to  100c,  but  the  experiment 
was  in  a  relatively  small  wind  tunnel,  which  likely  had  some  influence.  The  NASA-Langley 
turbulence  resource  only  provides  the  results  data  from  CFL3D  as  the  FUN3D  code  had 
almost  identical  results. 


Fig.  5.18:  Grid  system  for  flow  over  a  NACA  4412  airfoil. 


The  strand-Cartesian  solver  extends  the  Cartesian  grid  500  chords  in  all  direction  to 
reduce  the  effect  from  the  farfield  boundaries.  Riemann  invariant  boundary  conditions  are 
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imposed  at  the  farfield.  Portions  of  the  Cartesian  grid  system  are  shown  in  Figure  5.18.  A 
baseline  grid  consisting  of  10  levels  of  refinement  and  12,876  cells  is  used,  shown  in  Figure 
5.18(a).  Additionally,  a  wake-capturing  Cartesian  grid  consisting  of  15  levels  of  refinement 
and  411,908  cells  was  used,  shown  in  Figure  5.18(b).  The  wake  refining  grid  was  made  quite 
wide  to  handle  the  high  angle  of  attack.  The  total  number  of  cells  needed  could  be  reduced 
with  Cartesian  grid  adaption  strategies. 


Figure  5.19  shows  a  contour  plot  of  the  velocity  magnitude  over  the  airfoil.  Good  wake 
resolution  can  be  seen.  The  surface  pressure  coefficient  for  a  =  13.87°  is  shown  in  Figure 
5.20,  and  compared  with  the  experimental  data  of  Coles  and  Wadcock  [49],  and  the  CFD 
code  CFL3D.  Excellent  agreement  can  be  seen  between  the  strand-Cartesian  solver,  CFL3D 
and  the  experimental  data  of  Coles  and  Wadcock. 

Velocity  profiles  at  six  different  locations  along  the  trailing  edge  of  the  upper  surface 
are  shown  in  Figure  5.21.  Relatively  good  agreement  is  found  between  CFL3D  and  the 
strand-Cartesian  solver.  The  velocity  profiles  at  x  =  0.9528  and  x  =  0.8973  show  reverse 
flow  close  to  surface  of  the  airfoil  in  both  the  x  and  y  directions,  which  is  indicative  of 
the  flow  separation,  possibly  the  formation  of  a  laminar  separation  bubble.  More  accurate 
results  may  possibly  be  obtained  if  preconditioning  were  to  be  implemented.  The  trailing 
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x/c 


Fig.  5.20:  Pressure  coefficient  for  flow  over  a  NACA  4412  airfoil  at  M  =  0.09  and  Re  = 
1.52  x  106  at  an  angle  of  attack  of  13.87°. 


edge  of  the  airfoil  encounters  very  low  speed  due  to  the  separation  of  flow.  Preconditioning 
would  help  by  accelerating  the  convergence  to  a  steady  state  for  problems  where  a  significant 
portion  of  the  flow  is  low  speed. 

Table  5.4  shows  the  coefficient  of  lift  and  drag  for  CFL3D,  FUN3D  and  the  strand- 
Cartesian  solver.  The  strand-Cartesian  solver  predicted  a  coefficient  of  lift  close  to  FUN3D 
and  CFL3D.  The  drag  prediction  showed  some  discrepancy  with  FUN3D  and  CFL3D.  It  is 
possible  to  attribute  this  inaccuracy  to  the  lack  of  preconditioning  in  the  solver. 

Figure  5.22  shows  a  direction  field  plot  of  the  trailing  edge.  The  arrows  have  been  made 
a  uniform  size,  rather  than  sized  by  the  velocity  magnitude,  in  order  to  make  visualization 
easier.  Reverse  flow  in  the  plot  is  clearly  evident,  and  prediction  of  the  formation  of  a 
laminar  separation  bubble  on  the  trailing  edge  is  confirmed. 


Table  5.4:  Comparison  of  computed  lift  and  drag  coefficients  for  flow  over  a  NACA  4412 
airfoil,  at  M  =  0.09,  a  =  13.87°  Re  =  1.52  x  106. _ 


Ci 

cd 

Strand  (513  points) 

1.69013 

0.031477 

FUN3D  (513  points) 

1.7210 

0.02947 

CFL3D  (513  points) 

1.7170 

0.02861 
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Fig.  5.22:  Direction  field  plot  of  the  trailing  edge  separation  over  the  NACA  4412  airfoil 
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(b)  V  Velocity 

Fig.  5.21:  Velocity  profiles  along  the  NACA  4412  trailing  edge 


Chapter  6 
Conclusions 


6.1  Conclusions 

The  strand-Cartesian  grid  method  is  examined  for  its  ability  to  compute  high  Reynolds 
number  flows  over  a  variety  of  geometry  using  the  Spalart-Allmaras  (SA)  turbulence  model. 
Verification  of  the  model  is  undertaken  to  ensure  correct  implementation.  The  Method  of 
Manufactured  Solutions  grid  convergence  verification  method  showed  the  desired  second 
order  convergence  was  achieved  for  all  system  equations. 

The  validation  procedure  was  carried  out  by  comparing  results  from  the  strand-Cartesian 
solver  to  two  independent  compressible  codes;  FUN3D  and  CFL3D.  The  zero  pressure  gra¬ 
dient  flat  plate  was  run  at  Re  =  5  x  106.  The  strand  solver  computed  drag  to  within 
0.7%  and  1.5%  of  a  finer-grid  FUN3D  and  CFL3D  case  respectively.  Strand  grid  turbulent 
viscosity  and  velocity  profiles  at  various  locations  along  the  plate  were  shown  to  agree  well 
with  the  FUN3D  and  CFL3D. 

The  bump-in-channel  case  explored  Re  =  3  x  10b  flow  through  a  channel  with  a  pressure 
gradient  caused  by  a  bump.  The  strand  solver  computed  drag  to  within  0.3%  and  0.09% 
of  a  finer-grid  FUN3D  and  CFL3D  case  respectively.  Strand  grid  turbulent  viscosity  and 
velocity  profiles  at  various  locations  on  the  bump  were  shown  to  agree  excellently  with  the 
FUN3D  and  CFL3D.  Surface  pressure  and  skin  friction  coefficient  along  the  bump  was  also 
shown  to  agree  excellently  with  the  FUN3D  and  CFL3D. 

NACA  0012  case  examined  several  candidate  meshing  strategies,  including  strand 
smoothing  and  multi-strands,  with  and  without  telescoping  Cartesian  refinement  into  sharp 
corner  regions.  For  turbulent  Re  =  5  x  106  flow  over  a  NACA  0012  airfoil,  smooth  strands 
with  local  Cartesian  refinement  around  the  sharp  trailing  edge  computed  lift  to  within  0.1% 
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of  a  finer-grid  FUN3D  case,  and  drag  to  within  1.3%.  The  best  multi-strand  result  was  a 
difference  of  0.5%  in  lift  and  5.7%  in  drag. 

NACA  4412  case  examines  strand-Cartesian’s  ability  to  predict  flow  separation  on  the 
trailing  edge  of  an  airfoil.  For  turbulent  Re  =  1.52  x  106  flow  with  a  Mach  number  of  0.09 
over  a  NACA  4412  airfoil,  smooth  strands  with  local  Cartesian  refinement  around  the  sharp 
trailing  edge  computed  lift  to  within  1.79%  and  drag  to  within  6.8%  for  a  FUN3D  case, 
and  lift  to  within  1.56%  and  drag  to  within  10%  for  a  CFL3D  case.  Velocity  profiles  show 
reasonable  agreement  with  CFL3D,  however  implementing  preconditioning  to  the  solver  in 
the  future  may  increase  the  accuracy  of  the  solution. 

6.2  Future  Work 

The  Spalart-Allmaras  model  has  been  shown  to  excel  in  the  strand-Cartesian  solver 
with  the  given  cases  studies.  Further  cases  such  as  a  “backwards  facing  step”  should  be 
explored.  However,  cases  such  as  these  provide  a  challenge  in  meshing  with  strand  grids. 
Three-dimensional  cases  studies  should  be  performed,  however  the  challenges  with  strand 
meshing  in  3D  must  be  addressed  first.  Both  internal  and  external  sharp  corner  treatments 
must  be  examined  further.  For  example,  a  “3D  Supersonic  Square  Duct”  validation  case 
requires  a  highly  resolved  boundary  layer  on  adjacent  walls  of  a  square  duct,  causing  mesh 
generation  issues  at  the  internal  corner.  Such  challenges  need  be  overcome  before  a  strand- 
Cartesian  solver  may  be  validated  for  such  a  case. 

The  SA  model,  while  excellent  with  external  flows,  has  been  found  to  perform  poorly 
(or  at  least  worse  than  other  models)  with  interior  flows;  such  as  pipe  flow  [20] .  It  would  be 
appropriate  then,  to  implement  a  turbulence  model  that  excels  in  scenarios  that  SA  model 
struggles  with.  One  such  model  is  the  k  —  to  Menter  SST  model.  Like  the  SA  model,  this 
is  a  common  turbulence  model  with  plenty  of  documentation,  and  its  own  advantages  and 
disadvantages.  Unlike  the  SA  model,  the  k  —  oj  Menter  SST  model  has  been  shown  to  excel 
in  interior  flows  [51].  Verification  and  validation  should  be  undertaken  in  similar  fashion  to 
this  work. 
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Alternatively,  Large-Eddy  Simulation  (LES)  could  be  implemented  into  the  strand- 
Cartesian  solver  to  provide  highly  accurate  turbulence  simulations.  An  LES  strand-Cartesian 
solver  could  prove  to  be  extremely  powerful.  Other  work  that  could  be  explored  includes 
developing  high-order  methods  for  strands,  or  implementing  automated  grid  adaptation 
methods. 
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